Workflow:
Explicação de origem e relevância dos dados
Preparação e pré-processamento dos dados
Sumarização dos dados e gráficos
Análise estatística univariada
Análise de expressão diferencial e enriquecimento
1. Explicação de origem e relevância dos dados
Origem dos dados
Os dados provêm do estudo “LGG TCGA PanCancer Atlas 2018”, disponível no portal cBioPortal: 🔗 cBioPortal.
Esse estudo faz parte do The Cancer Genome Atlas (TCGA), um grande projeto colaborativo que visa caracterizar geneticamente diversos tipos de cancro usando tecnologias de alto rendimento, como RNA-Seq para expressão genética.
Tipo de cancro estudado
Este dataset está focado no Glioma de baixo grau (Low Grade Glioma - LGG), um tipo de tumor cerebral com evolução mais lenta do que o glioblastoma, mas que pode ser fatal em vários casos.
Apesar da sua progressão mais lenta, o LGG apresenta uma elevada heterogeneidade molecular, tornando-se relevante para estudos de estratificação de pacientes e identificação de subtipos tumorais.
Dados disponíveis
No cBioPortal, este estudo inclui:
- Dados clínicos (idade, sexo, sobrevida, etc.)
- Dados de Expressão Génica (RNA-Seq)
- Dados de Mutação Somática
- Dados de CNV (alterações no número de cópias)
- Dados de Metilação
- Dados de Fusões
Todos esses dados são *integráveis, permitindo análises multi-ómicas para uma visão mais ampla dos mecanismos moleculares envolvidos no desenvolvimento e progressão do **glioma*.
Dados de Expressão Génica (RNA-Seq)
Os dados de expressão génica foram processados usando o método *RSEM, com **normalização por lote* (batch-normalized) a partir da plataforma Illumina HiSeq RNASeqV2.
A matriz de expressão contém milhares de genes (geralmente
cerca de 20.000), permitindo análises em larga escala de:
- Expressão diferencial
- Coexpressão
- Redes regulatórias
Relevância científica
Este tipo de dados permite identificar genes diferencialmente
expressos entre:
- Subtipos tumorais
- Grupos com ou sem mutações específicas
- Amostras normais e tumorais
No caso do estudo, estão disponíveis apenas as 514 amostras tumorais de glioma de baixo grau (LGG).
Essas análises podem contribuir para:
- Descoberta de biomarcadores moleculares
- Estratificação de pacientes com base em perfis de
expressão
- Identificação de possíveis alvos terapêuticos
Além disso, os dados podem ser utilizados em análises de enriquecimento funcional (como GSEA ou over-representation analysis) para responder a perguntas como:
Quais as vias biológicas que estão mais ativas ou reprimidas em determinados grupos de pacientes?
Também é possível associar perfis de expressão a desfechos clínicos, como tempo de sobrevida ou resposta ao tratamento, contribuindo para avanços na medicina personalizada.
O dataset é compatível com técnicas de machine learning,
permitindo:
- Classificação de pacientes
- Seleção de features relevantes
- Predição de prognóstico
O acesso aberto e a documentação clara disponíveis pelo cBioPortal reforçam a reprodutibilidade científica e *a integração com ferramentas computacionais, como **APIs e bibliotecas em R ou Python*.
2. Preparação e pré-processamento dos dados
Workflow Dados Clínicos:
Visualização da estrutura dos dados
Remoção de NAs
Manipulação de dados
Sumarização e visualização de dados
Workflow Dados Expressão Génica:
Visualização estrutura dos dados
Remoção de duplicados
Tratamento de dados
Remoção de genes pouco expressos
Filtragem por níveis de expressão
Sumarização e visualização de dados
2.0. Importação das packages e dos dados
cran_packages <- c(
"ggplot2", # Visualization
"gplots", # Enhanced plots
"gridExtra", # Combine multiple plots
"DT", # Interactive datatables
"dplyr", # Data manipulation
"plyr", # Data manipulation (older)
"stringr", # String operations
"knitr", # RMarkdown rendering
"data.table", # Efficient data handling
"viridis", # Color palettes
"openxlsx", # Excel export
"corrplot", # Correlation plots
"mgcv", # GAM models
"ggpubr", # Publication-ready plots
"GGally", # Pair plots (extension of ggplot2)
"htmltools", # Plots
"caret", # Modelation
"randomForest", # Random Forest
"smotefamily", # SMOTE
"xgboost", # xgboost
"e1071", # SVM
"pROC", # ROC CURVE
"MLmetrics", # F1_score
"pander", # Tables
"tibble", # Interactive datatables
"Rtsne", # UMAP
"uwot", # T-SNE
"factoextra" # T-SNE
)
bioc_packages <- c(
"edgeR", # RNA-seq analysis
"limma", # Linear models for microarray/RNA-seq
"clusterProfiler", # Functional enrichment
"org.Hs.eg.db", # Human gene annotation
"AnnotationDbi", # Annotation infrastructure
"biomaRt", # Interface to BioMart databases
"enrichplot" # Visualization of enrichment results
)
other_packages <- c(
"gprofiler2" # g:Profiler interface (CRAN but bio-focused)
)
# Installation Block
# Install BiocManager if missing
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
# Install missing CRAN packages
cran_missing <- cran_packages[!cran_packages %in% rownames(installed.packages())]
if (length(cran_missing) > 0) install.packages(cran_missing)
# Install missing Bioconductor packages
bioc_missing <- bioc_packages[!bioc_packages %in% rownames(installed.packages())]
if (length(bioc_missing) > 0) BiocManager::install(bioc_missing)
# Install other packages (gprofiler2 is from CRAN)
other_missing <- other_packages[!other_packages %in% rownames(installed.packages())]
if (length(other_missing) > 0) install.packages(other_missing)
# Load All Packages
all_packages <- c(cran_packages, bioc_packages, other_packages)
invisible(lapply(all_packages, function(pkg) {
suppressMessages(library(pkg, character.only = TRUE))
}))
# Import data
clini_data <- read.delim("lgg_tcga_pan_can_atlas_2018_clinical_data.tsv", stringsAsFactors=TRUE)
gene_data <- read.delim("all_genes_mrna.txt")Inicialmente, foram importados os dados de expressão génica,
gene_data, e os dados clínicos dos pacientes analisados,
clini_data.
2.1. Dados Clínicos
2.1.1. Estrutura dos dados
A estrutura dos dados clínicos é apresentada em baixo. Os pacientes
são identificados por um ID (Patient.ID) que faz
correspondência com os dados de expressão génica. Nesta fase estão
presentes 63 variáveis e 514 observações (listado abaixo).
## 'data.frame': 514 obs. of 63 variables:
## $ Study.ID : Factor w/ 1 level "lgg_tcga_pan_can_atlas_2018": 1 1 1 1 1 1 1 1 1 1 ...
## $ Patient.ID : Factor w/ 514 levels "TCGA-CS-4938",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Sample.ID : Factor w/ 514 levels "TCGA-CS-4938-01",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Diagnosis.Age : int 31 67 44 37 50 47 39 40 43 53 ...
## $ Neoplasm.Disease.Stage.American.Joint.Committee.on.Cancer.Code : logi NA NA NA NA NA NA ...
## $ American.Joint.Committee.on.Cancer.Publication.Version.Type : logi NA NA NA NA NA NA ...
## $ Aneuploidy.Score : int 1 7 2 5 1 2 0 1 8 6 ...
## $ Buffa.Hypoxia.Score : int -27 -33 -27 -11 -11 -27 -25 -31 -21 -35 ...
## $ Cancer.Type : Factor w/ 1 level "Glioma": 1 1 1 1 1 1 1 1 1 1 ...
## $ TCGA.PanCanAtlas.Cancer.Type.Acronym : Factor w/ 1 level "LGG": 1 1 1 1 1 1 1 1 1 1 ...
## $ Cancer.Type.Detailed : Factor w/ 4 levels "Astrocytoma",..: 1 1 1 1 1 4 1 1 4 4 ...
## $ Last.Communication.Contact.from.Initial.Pathologic.Diagnosis.Date : int 3574 NA NA 552 1828 1966 1222 8 NA 1631 ...
## $ Birth.from.Initial.Pathologic.Diagnosis.Date : int -11509 -24578 -16297 -13565 -18494 -17460 -14418 -14920 -15848 -19399 ...
## $ Last.Alive.Less.Initial.Pathologic.Diagnosis.Date.Calculated.Day.Value : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Disease.Free..Months. : num NA NA NA NA NA ...
## $ Disease.Free.Status : Factor w/ 2 levels "0:DiseaseFree",..: NA NA NA NA NA 1 1 NA NA NA ...
## $ Months.of.disease.specific.survival : num 117.5 7.69 43.89 36.36 60.1 ...
## $ Disease.specific.Survival.status : Factor w/ 2 levels "0:ALIVE OR DEAD TUMOR FREE",..: 1 2 2 2 1 1 1 1 2 1 ...
## $ Ethnicity.Category : Factor w/ 2 levels "Hispanic Or Latino",..: 2 2 NA NA NA NA NA NA NA 2 ...
## $ Form.completion.date : Factor w/ 170 levels "1/10/12","1/13/14",..: 63 97 83 84 98 84 84 95 82 85 ...
## $ Fraction.Genome.Altered : num 0.0518 0.2241 0.0937 0.1625 0.0603 ...
## $ Genetic.Ancestry.Label : Factor w/ 8 levels "AFR","AFR_ADMIX",..: 5 5 1 5 5 5 5 5 1 5 ...
## $ Neoplasm.Histologic.Grade : Factor w/ 2 levels "G2","G3": 1 2 2 2 1 1 2 2 1 2 ...
## $ Neoadjuvant.Therapy.Type.Administered.Prior.To.Resection.Text : Factor w/ 3 levels "No","Yes","Yes (Pharmaceutical Treatment Prior To Resection)": 1 1 1 1 1 1 1 1 1 1 ...
## $ ICD.10.Classification : Factor w/ 6 levels "C71.0","C71.1",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ International.Classification.of.Diseases.for.Oncology..Third.Edition.ICD.O.3.Histology.Code: Factor w/ 5 levels "9382/3","9400/3",..: 2 3 3 3 2 4 3 3 1 5 ...
## $ International.Classification.of.Diseases.for.Oncology..Third.Edition.ICD.O.3.Site.Code : Factor w/ 6 levels "C71.0","C71.1",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ In.PanCan.Pathway.Analysis : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ Informed.consent.verified : Factor w/ 1 level "Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ MSI.MANTIS.Score : num 0.303 0.274 0.281 0.275 0.27 ...
## $ MSIsensor.Score : num 0 0 0.02 0.25 0.04 0.1 0 0.16 0.07 0 ...
## $ Mutation.Count : int 14 43 25 24 21 43 24 22 42 28 ...
## $ New.Neoplasm.Event.Post.Initial.Therapy.Indicator : Factor w/ 2 levels "No","Yes": NA 2 2 NA 1 NA 1 NA 2 NA ...
## $ Oncotree.Code : Factor w/ 4 levels "DIFG","LGGNOS",..: 1 1 1 1 1 4 1 1 4 4 ...
## $ Overall.Survival..Months. : num 117.5 7.69 43.89 36.36 60.1 ...
## $ Overall.Survival.Status : Factor w/ 2 levels "0:LIVING","1:DECEASED": 1 2 2 2 1 1 1 1 2 1 ...
## $ Other.Patient.ID : Factor w/ 513 levels "001AD307-4AD3-4F1D-B2FC-EFC032871C7E",..: 98 501 77 295 205 401 471 301 493 365 ...
## $ American.Joint.Committee.on.Cancer.Metastasis.Stage.Code : logi NA NA NA NA NA NA ...
## $ Neoplasm.Disease.Lymph.Node.Stage.American.Joint.Committee.on.Cancer.Code : logi NA NA NA NA NA NA ...
## $ American.Joint.Committee.on.Cancer.Tumor.Stage.Code : logi NA NA NA NA NA NA ...
## $ Person.Neoplasm.Cancer.Status : Factor w/ 2 levels "Tumor Free","With Tumor": NA 2 2 2 2 1 1 NA 2 2 ...
## $ Progress.Free.Survival..Months. : num 117.5 0.296 38.926 36.361 60.098 ...
## $ Progression.Free.Status : Factor w/ 2 levels "0:CENSORED","1:PROGRESSION": 1 2 2 2 1 1 1 1 2 1 ...
## $ Primary.Lymph.Node.Presentation.Assessment : logi NA NA NA NA NA NA ...
## $ Prior.Diagnosis : Factor w/ 3 levels "No","Yes","Yes, History Of Synchronous And Or Bilateral Malignancy": 1 1 1 1 1 1 1 1 1 1 ...
## $ Race.Category : Factor w/ 4 levels "American Indian or Alaska Native",..: 4 4 3 4 4 4 4 4 3 4 ...
## $ Radiation.Therapy : Factor w/ 2 levels "No","Yes": 1 2 2 1 2 2 2 1 2 2 ...
## $ Ragnum.Hypoxia.Score : int -22 -22 -22 6 -16 -12 -20 -24 -18 -20 ...
## $ Number.of.Samples.Per.Patient : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Sample.Type : Factor w/ 1 level "Primary": 1 1 1 1 1 1 1 1 1 1 ...
## $ Sex : Factor w/ 2 levels "Female","Male": 1 2 1 2 2 1 2 2 2 1 ...
## $ Somatic.Status : Factor w/ 1 level "Matched": 1 1 1 1 1 1 1 1 1 1 ...
## $ Subtype : Factor w/ 3 levels "LGG_IDHmut-codel",..: 2 3 2 2 2 1 2 2 3 1 ...
## $ Tumor.Break.Load : int 31 26 22 130 20 6 49 5 28 14 ...
## $ Tissue.Prospective.Collection.Indicator : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Tissue.Retrospective.Collection.Indicator : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ Tissue.Source.Site : Factor w/ 26 levels "Asterand","Case Western",..: 21 21 21 21 21 21 21 21 21 21 ...
## $ Tissue.Source.Site.Code : Factor w/ 26 levels "CS","DB","DH",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ TMB..nonsynonymous. : num 0.467 1.467 0.867 0.8 0.7 ...
## $ Tumor.Disease.Anatomic.Site : Factor w/ 1 level "CNS": 1 1 1 1 1 1 1 1 1 1 ...
## $ Tumor.Type : Factor w/ 4 levels "Astrocytoma",..: 1 1 1 1 1 4 1 1 4 4 ...
## $ Patient.Weight : logi NA NA NA NA NA NA ...
## $ Winter.Hypoxia.Score : int -38 -46 -32 -26 -22 -30 -32 -38 -20 -44 ...
2.1.2. Manipulação de dados
A manipulação dos dados clínicos consistiu nos seguintes processos:
- Conversão do `Patient.ID` para o mesmo formato do dataset
gene_data; - Remoção de variáveis clínicas consideradas menos relevantes no âmbito do estudo - este processo foi validado com literatura associada aos dados;
- Simplificação dos nomes das colunas;
- Conversão variável de meses de sobrevivência para anos de
sobrevivência -
surv_years;
# Change column names
nomes <- names(clini_data)
nomes <- gsub("[\\.]", " ", nomes)
nomes <- gsub(" ", " ", nomes)
nomes <- gsub(" $", "", nomes, perl=T)
names(clini_data) <- nomes
# Filter by columns of interest
novas_colunas<-c("Sample ID","Diagnosis Age","Cancer Type Detailed","Months of disease specific survival","Disease specific Survival status","Fraction Genome Altered","Genetic Ancestry Label","Sex","Tumor Break Load")
clini_data <- clini_data[,novas_colunas]
# Change variable names
colnames(clini_data) <- c("sample_ID","diagnosis_age","cancer_type","surv_months","surv_status","genome_alt","ancestry","sex","tbl")
# Convert sample ID to same format as expression data
sample_id <- clini_data$sample_ID
sample_id <- gsub("-", "\\.", sample_id)
clini_data$sample_ID <- sample_id
# Convert survival months to years
clini_data$surv_months <- clini_data$surv_months/12
colnames(clini_data) <- c("sample_ID","diagnosis_age","cancer_type","surv_years","surv_status","genome_alt","ancestry","sex","tbl")
2.1.3. Limpeza de Dados - Remoção de NAs
Uma vez que a proporção de NAs é baixa (6%), optámos por remover as linhas com NAs (Listwise Deletion). Por consequente, a amostra ficou reduzida a 481 pacientes.
# Remoção de NAs
clini_data_no_na <- na.omit(clini_data)
# Proporção de NAs no dataset
(nrow(clini_data) - nrow(clini_data_no_na))/nrow(clini_data)## [1] 0.06420233
2.1.4. Sumarização e Visualização de Dados
## sample_ID diagnosis_age cancer_type
## Length:481 Min. :14.00 Astrocytoma :183
## Class :character 1st Qu.:33.00 Low-Grade Glioma (NOS): 0
## Mode :character Median :41.00 Oligoastrocytoma :125
## Mean :43.02 Oligodendroglioma :173
## 3rd Qu.:53.00
## Max. :87.00
##
## surv_years surv_status genome_alt
## Min. : 0.000 0:ALIVE OR DEAD TUMOR FREE:371 Min. :0.0000
## 1st Qu.: 1.090 1:DEAD WITH TUMOR :110 1st Qu.:0.0517
## Median : 1.797 Median :0.0891
## Mean : 2.638 Mean :0.1145
## 3rd Qu.: 3.348 3rd Qu.:0.1503
## Max. :17.597 Max. :0.9477
##
## ancestry sex tbl
## EUR :436 Female:213 Min. : 0.00
## AFR : 15 Male :268 1st Qu.: 5.00
## AFR_ADMIX: 8 Median : 16.00
## EAS : 8 Mean : 32.66
## EUR_ADMIX: 7 3rd Qu.: 37.00
## AMR : 4 Max. :618.00
## (Other) : 3O resumo estatístico dos dados permite conferir que não existem anomalias nos dados (e.g. proporções acima de 100% ou número de anos negativos). Seguem abaixo observações sobre as variáveis:
Diagnosis age - vasta distribuição, variando entre 14 e 87 anos;
Cancer type - Oligoastrocytoma representa 25% dos dados, sendo que os restantes são mais frequentes - Oligodendroglioma (36%) e Astrocytoma (38%);
Survival years - 75% dos dados variam entre 0 e 3.4 anos, com um máximo de 17.6 anos, indicando que a variável poderá conter outliers
Survival status - 70% da amostra apresenta-se em remissão (vivo ou não);
Genome altered - a fração de genoma alterado apresenta uma grande variação, entre 0% e 95% de alteração. No entanto, sendo que 75% dos dados não ultrapassam 15% de alteração, valores muito elevados poderão ser outliers;
Ancestry - a maioria da amostra (91%) tem ancestralidade europeia;
Sex - a amostra encontra-se sensivelmente equilibrada, com 44% de mulheres e 56% de homens;
Tumor Break Load (TBL) - esta métrica avalia a instabilidade genómica; 75% dos dados variam entre 0 e 37, enquanto que o valor máximo é de 618. À semelhança de outras variáveis, esta também aparenta ter outliers.
Numerical variables:
Diagnosis Age (
diagnosis_age)Years of survival (
surv_years)Fraction Genome Altered (
genome_alt)Tumor Break Load (
tbl)
Factor variables:
Sample ID (
sample_ID)Sex (
sex)Cancer Type (
cancer_type)Genetic Ancestry (
ancestry)Survival Status (
surv_status)
Visualização univariada
# Variáveis numéricas
a <- ggplot(alpha=0.8) + geom_histogram(data=clini_data_no_na, aes(x=diagnosis_age), fill="#7eb0d5") + theme_bw()
b <- ggplot(alpha=0.8) + geom_histogram(data=clini_data_no_na, aes(x=surv_years), fill="#7eb0d5") + theme_bw()
ggarrange(a, b)
A variável Diagnosis Age tem uma distribuição próxima da distribuição normal, enquanto que a Survival Years tem uma distribuição enviesada à direita.
# Variáveis numéricas
a <- ggplot(alpha=0.8) + geom_histogram(data=clini_data_no_na, aes(x=genome_alt), fill="#7eb0d5") + theme_bw()
b <- ggplot(alpha=0.8) + geom_histogram(data=clini_data_no_na, aes(x=tbl), fill="#7eb0d5") + theme_bw()
grid.arrange(a, b, ncol=2)
Ambas as variáveis apresentam uma distribuição enviesada à direita.
# Variáveis categóricas
a <- ggplot(alpha=0.8) + geom_bar(data=clini_data_no_na, aes(x=sex, fill=sex)) + scale_fill_manual(values = palette_3) + theme_bw()
b <- ggplot(alpha=0.8) + geom_bar(data=clini_data_no_na, aes(x=cancer_type, fill=cancer_type)) + scale_fill_manual(values = palette_3) + theme_bw()
grid.arrange(a, b, ncol=2)
A variável categóriga Sex está distribuida de forma equilibrada nos dois géneros. A variável Cancer Type apresenta uma pequena variação entre os 3 tipos de cancro como referido anteriormente.
# Variáveis categóricas
a <- ggplot(alpha=0.8) + geom_bar(data=clini_data_no_na, aes(x=ancestry, fill=ancestry)) + scale_fill_manual(values = palette_3) + theme_bw()
b <- ggplot(alpha=0.8) + geom_bar(data=clini_data_no_na, aes(x=surv_status, fill=surv_status)) + scale_fill_manual(values = palette_3) + theme_bw()
grid.arrange(a, b, ncol=2)
A análise da variável Ancestry confirma o que foi dito anteriormente, apresentando a ancestralidade europeia como a amostra predominante. A variável Survival Status demonstra uma taxa elevada de pacientes em remissão (vivos ou não) em comparação com os que morreram com tumor.
Visualização multivariada
Variável numérica vs. numérica
Os pares Diagnosis Age - Survival Years e Diagnosis Age - Genome
Alteration têm uma correlação baixa, negativa e positiva,
respetivamente.
Os pares Diagnosis Age - Tumor Break Load, Survival Years - Genome
Alteration and Survival Years - Tumor Break Load apresentam uma linha
praticamente horizontal, pelo que não aparentam ter correlação.
O par Genome Alteration - Tumor Break Load apresenta correlação na ordem
dos 0.484, que se mantém abaixo do coeficiente de correlação de 0.5,
pelo que não aparenta haver multicolinearidade.
Variável numérica vs. categórica
# Remover outliers
# get outliers
out_age <- which(clini_data_no_na$diagnosis_age>(quantile(clini_data_no_na$diagnosis_age, .75) + 1.5*IQR(clini_data_no_na$diagnosis_age)))
out_surv <- which(clini_data_no_na$surv_years>(quantile(clini_data_no_na$surv_years, .75) + 1.5*IQR(clini_data_no_na$surv_years)))
out_gen <- which(clini_data_no_na$genome_alt>(quantile(clini_data_no_na$genome_alt, .75) + 1.5*IQR(clini_data_no_na$genome_alt)))
out_tbl <- which(clini_data_no_na$tbl>(quantile(clini_data_no_na$tbl, .75) + 1.5*IQR(clini_data_no_na$tbl)))
# remove outliers
no_out_data <- clini_data_no_na[-c(out_age, out_surv, out_gen, out_tbl),]Após uma análise grafica inicial constatou-se que algumas variáveis tinham outliers, pelo que se realizou a remoção dos mesmos.
Variável sex
#### sex vs. diagnosis age
a <- ggplot(alpha=0.8, data=no_out_data, aes(x=sex, y=diagnosis_age, fill=sex)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
#### sex vs. genome altered
b <- ggplot(alpha=0.8, data=no_out_data, aes(x=sex, y=genome_alt, fill=sex))+ geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
ggarrange(a, b, common.legend=T)
A distribuição das variáveis Diagnosis Ages e Genome Alteration em
função da variável Sex, representadas acima, não apresentam grande
variação entre sexos. Infere-se que não há uma associação entre o sexo e
essas variáveis.
#### sex vs. tbl
a <- ggplot(alpha=0.8, data=no_out_data, aes(x=sex, y=tbl, fill=sex)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
#### sex vs. survival
b <- ggplot(alpha=0.8, data=no_out_data, aes(x=sex, y=surv_years, fill=sex)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
ggarrange(a, b, common.legend=T)
A distribuição das variáveis Tumor Break Load e Survival Years em função da variável Sex, representadas acima, não apresentam grande variação entre sexos. Infere-se que não há uma associação entre o sexo e essas variáveis.
Variável Cancer Type
#### cancer type vs. diagnosis age
ggplot(alpha=0.8, data=no_out_data, aes(x=cancer_type, y=diagnosis_age, fill=sex)) + geom_violin() + geom_boxplot(width=.1,position = position_dodge(width =0.9)) + scale_fill_manual(values = palette_3) + theme_bw()
Neste gráfico, o cancro Oligodendroglioma apresenta uma distribuição com mediana dos valores de Diagnosis Age acima dos outros tipos de cancro.
#### cancer type vs. genome altered
ggplot(alpha=0.8, data=no_out_data, aes(x=cancer_type, y=genome_alt, fill=sex)) + geom_violin() + geom_boxplot(width=.1,position = position_dodge(width =0.9)) + scale_fill_manual(values = palette_3) + theme_bw()
Neste gráfico, foi feita uma distinção entre sexos, devido à ligeira
diferença observada no gráfico individual de Genome Altered. Observa-se
pouca variação entre os 3 tipos de cancro em relação à alteração de
genoma. Entre sexos dentro dos pares há uma variação mais notável nos
pacientes com oligoastrocytoma.
#### cancer type vs. tbl
a <- ggplot(alpha=0.8, data=no_out_data, aes(x=cancer_type, y=tbl, fill=cancer_type)) + geom_violin() + geom_boxplot(width=.1,position = position_dodge(width =0.9)) + scale_fill_manual(values = palette_3) + theme_bw()
#### cancer type vs. years survival
b <- ggplot(alpha=0.8, data=no_out_data, aes(x=cancer_type, y=surv_years, fill=cancer_type)) + geom_violin() + geom_boxplot(width=.2,position = position_dodge(width =0.9)) + scale_fill_manual(values = palette_3) + theme_bw()
c <- ggplot(alpha=0.8, data=no_out_data, aes(x=cancer_type, y=surv_years, fill=cancer_type)) + geom_violin() + geom_boxplot(width=.2,position = position_dodge(width =0.9)) + scale_fill_manual(values = palette_3) + theme_bw()
ggarrange(a, b, common.legend=T)
No primeiro gráfico verifica-se uma grande variação dos valores de
Tumor Break Load entre tipos de cancro.
No segundo gráfico não se observa grande variação entre tipos de cancro
relativamente aos valores de Survival Years.
Variável Ancestry
#### ancestry vs. genome altered
a <- ggplot(no_out_data, aes(x = ancestry, y = genome_alt)) +
geom_boxplot(fill = palette_3[1:length(unique(no_out_data$ancestry))], alpha = 0.5) +
geom_jitter(width = 0.2, alpha = 0.8, color = "grey10") +
xlab("") +
theme_bw()
#### ancestry vs. age
b <- ggplot(no_out_data, aes(x = ancestry, y = diagnosis_age)) +
geom_boxplot(fill = palette_3[1:length(unique(no_out_data$ancestry))], alpha = 0.5) +
geom_jitter(width = 0.2, alpha = 0.8, color = "grey10") +
xlab("") +
theme_bw()
ggarrange(a, b, common.legend=T)
Em ambos os gráficos mostram diferenças significativas entre essas variáveis e essa variação pode ser devido a dimensões de amostra distintas.Embora dentro dos grupos com maiores dimensões (AFR e EUR) não há diferenças muito notáveis.
#### ancestry vs. tbl
a <- ggplot(no_out_data, aes(x = ancestry, y = tbl)) +
geom_boxplot(fill = palette_3[1:length(unique(no_out_data$ancestry))], alpha = 0.5) +
geom_jitter(width = 0.2, alpha = 0.8, color = "grey10") +
xlab("") +
theme_bw()
#### ancestry vs. survival
b <- ggplot(no_out_data, aes(x = ancestry, y = surv_years)) +
geom_boxplot(fill = palette_3[1:length(unique(no_out_data$ancestry))], alpha = 0.5) +
geom_jitter(width = 0.2, alpha = 0.8, color = "grey10") +
xlab("") +
theme_bw()
ggarrange(a, b, common.legend=T)
Como observado no caso acima, ambos os gráficos mostram diferenças significativas entre essas variáveis e essa variação pode ser devido a dimensões de amostra distintas.Embora dentro dos grupos com maiores dimensões (AFR e EUR) não há diferenças muito notáveis.
Variável surv_status
#### sex vs. diagnosis age
a <- ggplot(alpha=0.8, data=no_out_data, aes(x=surv_status, y=diagnosis_age, fill=surv_status)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
#### sex vs. genome altered
b <- ggplot(alpha=0.8, data=no_out_data, aes(x=surv_status, y=genome_alt, fill=surv_status))+ geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
ggarrange(a, b, common.legend=T)
No primeiro plot observa-se uma clara distinção entre os fatores de Survival Status com base na Diagnosis Age. Diagnosis Age tem uma mediana mais alta para Dead With Tumor. O mesmo se verifica entre a variável Genome Alteration e Survival Status.
#### sex vs. tbl
a <- ggplot(alpha=0.8, data=no_out_data, aes(x=surv_status, y=tbl, fill=surv_status)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
#### sex vs. survival
b <- ggplot(alpha=0.8, data=no_out_data, aes(x=surv_status, y=surv_years, fill=surv_status)) + geom_violin() + geom_boxplot(width=.1) + scale_fill_manual(values = palette_3) + theme_bw()
ggarrange(a, b, common.legend=T)
No primeiro plot observa-se uma clara distinção entre os fatores de Survival Status com base no Tumor Break Load. TBL tem uma mediana mais alta para Dead With Tumor. Verifica-se uma variação menos notável entre a variável Survival Years e Survival Status.
2.2. Dados de Expressão Génica
2.2.1. Estrutura dos dados
A estrutura dos dados de expressão génica é apresentada em baixo. Os pacientes são identificados por um ID (Patient.ID) que faz correspondência com os dados clínicos. Nesta fase estão presentes 10 variáveis e 20531 observações (listado abaixo).
## Hugo_Symbol Entrez_Gene_Id TCGA.CS.4938.01 TCGA.CS.4941.01 TCGA.CS.4942.01
## 1 100130426 0.0000 0.0000 0.0000
## 2 100133144 8.7141 36.4493 11.8131
## 3 UBE2Q2P2 100134869 22.7523 21.1767 11.0242
## 4 HMGB1P1 10357 268.5760 156.6870 185.1380
## 5 10431 845.8150 390.2690 621.4530
## 6 136542 0.0000 0.0000 0.0000## 'data.frame': 20531 obs. of 10 variables:
## $ Hugo_Symbol : chr "" "" "UBE2Q2P2" "HMGB1P1" ...
## $ Entrez_Gene_Id : int 100130426 100133144 100134869 10357 10431 136542 155060 26823 280660 317712 ...
## $ TCGA.CS.4938.01: num 0 8.71 22.75 268.58 845.82 ...
## $ TCGA.CS.4941.01: num 0 36.4 21.2 156.7 390.3 ...
## $ TCGA.CS.4942.01: num 0 11.8 11 185.1 621.5 ...
## $ TCGA.CS.4943.01: num 0 8.61 5.08 269.84 835.73 ...
## $ TCGA.CS.4944.01: num 0 0 30.3 216.3 812.5 ...
## $ TCGA.CS.5390.01: num 0 5.34 27.89 159.76 576.9 ...
## $ TCGA.CS.5393.01: num 0 3.78 8.72 198.19 551.95 ...
## $ TCGA.CS.5394.01: num 0 8.31 15.45 208.54 607.9 ...
2.2.2. Remoção de duplicados e NAs
Nesta fase, procedeu-se a retirar genes duplicados, genes descontinuados e NAs do dataset. Em seguida, transpôs-se o dataset para colocar os Ids dos genes nas linhas, algo essencial para utilizar o package (EdgeR) nos passos seguintes. De modo a simplificar e ser possível comparar dados clínicos com dados de expressão génica excluiram-se os pacientes sem dados clínicos.
# Gene descontinuado DGCR9, row nº 4886 - será removido
gene_data_nd <- gene_data[-c(4886),]
# Remove duplicates
gene_data_nd <- gene_data_nd[-c(which(duplicated(gene_data_nd$Entrez_Gene_Id))),]
# Remoção de NAs
gene_data_nona <- na.omit(gene_data_nd)
# Proporção de NAs no dataset
(nrow(gene_data_nd) - nrow(gene_data_nona))/nrow(gene_data_nd)## [1] 0# Converter Gene ID em rownames
rownames(gene_data_nona) <- gene_data_nona$Entrez_Gene_Id
# Exclusão de pacientes sem dados clínicos
complete_ids <- clini_data_no_na$sample_ID
gene_data_nona <- gene_data_nona[,complete_ids]
2.2.3. Tratamento de dados
Para ser possível filtrar os dados usou-se DGEList e agrupou-se os dados em 4 grupos - f_cancer, f_survival, f_ancestry and f-sex.
# Converter em formato DGEList (package edgeR)
## sem agrupamento
d0 <- DGEList(gene_data_nona) # no grouping
# criação de grupos
f_cancer <- revalue(clini_data_no_na$cancer_type, c("Astrocytoma" = "A", "Oligoastrocytoma" = "B", "Oligodendroglioma" = "C", "Low-Grade Glioma (NOS)" = "D"))
f_survival <- revalue(clini_data_no_na$surv_status, c("0:ALIVE OR DEAD TUMOR FREE" = "A", "1:DEAD WITH TUMOR" = "B"))
f_ancestry <- revalue(clini_data_no_na$ancestry, c("EUR" = "A", "AFR" = "B", "AFR_ADMIX" = "C", "EAS" = "D", "EUR_ADMIX" = "E", "AMR" = "F", "SAS" = "G", "SAS_ADMIX" = "H"))
f_sex <- revalue(clini_data_no_na$sex, c("Male" = "A", "Female" = "B"))
# agrupamento por fatores
# d0_g <- DGEList(gene_data_nona, group=f_cancer) # or f_sex, etc.
# Visualizar estrutura dos dados
head(d0$samples)## group lib.size norm.factors
## TCGA.CS.4938.01 1 21134536 1
## TCGA.CS.4941.01 1 18815260 1
## TCGA.CS.4942.01 1 18777593 1
## TCGA.CS.4943.01 1 17861324 1
## TCGA.CS.4944.01 1 24990566 1
## TCGA.CS.5393.01 1 18182161 1
2.2.4. Remoção de genes pouco expressos
A primeira fase de filtragem passou por excluir os genes pouco expressos, reduzindo a lista de 20504 elementos para 14357.
## Warning in filterByExpr.DGEList(d0): All samples appear to belong to the same
## group.
2.2.5. Filtragem por níveis de expressão
Numa segunda fase de filtragem, utilizaram-se filtros flat patterns nomeadamente aplicar 2*mediana (Filtro 1) e max/min > 2 (Filtro 2), separadamente.
# Filtros flat patterns - 2*mediana
gene_exp_sd <- apply(gene_data_nona, 1, sd)
d_med <- 2*median(gene_exp_sd)
high_exp <- which(gene_exp_sd > d_med)
mean_exp <- data.frame(gene_id=high_exp,exp=apply(gene_data_nona[high_exp,], 1, mean), row.names = c())
mean_exp$exp_max <- mean_exp$exp/max(mean_exp$exp)
mean_exp$exp_log <- log(mean_exp$exp)# Filtros flat patterns - max/min > 2
gene_min <- apply(gene_data_nona, 1, min)
gene_max <- apply(gene_data_nona, 1, max)
rat <- which(gene_max/gene_min > 2)
mean_exp_f2 <- data.frame(id=rat, exp=apply(gene_data_nona[rat,-1], 1, mean))
mean_exp_f2$exp_max <- mean_exp_f2$exp/max(mean_exp_f2$exp)
mean_exp_f2$exp_log <- log(mean_exp_f2$exp)
2.2.6. Sumarização e visualização de dados
Filtro flat pattern 1
Recorrendo aos dados filtrados com o filtro 1, realizaram-se histogramas para visualizar a distribuição dos dados (Gráfico da esquerda). No gráfico da direita, apenas se fez uma normalização dos dados com log transformation de modo a poder visualizar melhor a distribuição, que aparenta ser próxima do normal.
# Filtro flat pattern 1
a <- ggplot(mean_exp) + geom_histogram(aes(x=exp), bins=400, fill="#7eb0d5") +ggtitle("Flat pattern 1")
b <- ggplot(mean_exp) + geom_histogram(aes(x=exp_log), bins=400, fill="#7eb0d5") + ggtitle("Flat pattern 1 - log transformation")
grid.arrange(a, b, ncol=2)
Filtro flat pattern 2
Recorrendo aos dados filtrados com o filtro 2, realizaram-se histogramas para visualizar a distribuição dos dados (Gráfico da esquerda). No gráfico da direita, apenas se fez uma normalização dos dados com log transformation de modo a poder visualizar melhor a distribuição, que aparenta estar enviesada à esquerda.
# Filtro flat pattern 2
a <- ggplot(mean_exp_f2) + geom_histogram(aes(x=exp), bins=200, fill="#7eb0d5") + ggtitle("Flat pattern 2")
b <- ggplot(mean_exp_f2) + geom_histogram(aes(x=exp_log), bins=200, fill="#7eb0d5") + ggtitle("Flat pattern 2 - log transformation")
grid.arrange(a, b, ncol=2)
Ontologia de genes
# 20 genes mais expressos
all_genes_means <- data.frame(gene_id=c(rownames(gene_data_nona)),exp=apply(gene_data_nona, 1, mean))
all_genes_means <- all_genes_means[order(all_genes_means$exp, decreasing=T),]
high_20 <- head(all_genes_means, n=20)
# Visualização das funções
go_prof <- gost(high_20$gene_id, organism = "hsapiens")
gostplot(go_prof, capped=F, interactive = T)O gráfico seguinte representa a ontologia dos 20 genes mais
expressos. Cada círculo representa uma correspondência de um gene a uma
das suas funções. Genes podem estar associados a mais do que uma função.
Duas das ontologias com maior representação pelos genes mais expressos
são de componentes celulares (GO:CC) e
processos biológicos (GO:BP).
3. Análise Univariada
# Preparar dados clínicos e fatores
clean_data_for_2part <- clini_data_no_na
f_ancestralidade <- as.factor(clean_data_for_2part$sex)
names(f_sex) <- clean_data_for_2part$sample_ID
f_cancer <- as.factor(clean_data_for_2part$cancer_type)
names(f_cancer) <- clean_data_for_2part$sample_ID
f_ancestralidade <- as.factor(clean_data_for_2part$ancestry)
names(f_ancestralidade) <- clean_data_for_2part$sample_ID
f_survival <- as.factor(clean_data_for_2part$surv_status)
names(f_survival) <- clean_data_for_2part$sample_ID
fatores <- list(
Sexo = f_sex,
TipoTumor = f_cancer,
Ancestralidade = f_ancestralidade,
Sobrevivencia = f_survival
)
# Expressão log2-CPM
log_cpm <- cpm(d0_f, log = TRUE)
# Inicializar variáveis
resultados_df <- data.frame()
plot_list <- list()
i <- 1
# Loop por genes e variáveis clínicas
for (gene_id in high_20$gene_id) {
for (nome_var in names(fatores)) {
grupo <- fatores[[nome_var]]
amostras_comuns <- intersect(colnames(log_cpm), names(grupo))
if (length(amostras_comuns) < 2) next
grupo_valid <- grupo[amostras_comuns]
expr <- log_cpm[as.character(gene_id), amostras_comuns]
if (length(unique(grupo_valid)) < 2) next
df_plot <- data.frame(expr = expr, grupo = grupo_valid)
# Boxplot
p <- ggplot(df_plot, aes(x = grupo, y = expr, fill=grupo)) +
geom_boxplot() +
scale_fill_manual(values = palette_3) + theme_bw() +
labs(title = paste(gene_id, "-", nome_var), x = nome_var, y = "Log2 CPM")
plot_list[[i]] <- p
i <- i + 1
# Teste estatístico
p_val <- if (length(unique(grupo_valid)) == 2) {
t.test(expr ~ grupo_valid)$p.value
} else {
summary(aov(expr ~ grupo_valid))[[1]][["Pr(>F)"]][1]
}
resultados_df <- rbind(resultados_df, data.frame(
gene = gene_id,
variavel_clinica = nome_var,
p_value = p_val
))
}
}
# Mostrar os boxplots em grupos de 4 com espaços para comentários
text_chunks <- c(
"Expressão muito elevada e bastante estável entre os grupos. Não há variações visíveis relevantes por sexo, tipo tumoral ou sobrevivência.",
"Apresenta expressão mais elevada em um subtipo tumoral específico. Discreta variação por ancestralidade.",
"Tipo tumoral mostra-se como principal fator diferenciador; as restantes variáveis não apresentam padrões marcantes.",
"Diferença clara entre subtipos tumorais. Leves variações por sexo e ancestralidade, mas não significativas.",
"Padrão consistente com maior expressão em um subtipo tumoral. Outras variáveis sem impacto aparente.",
"Variação moderada por tipo tumoral e ancestralidade. Não se observa tendência clara por sexo ou sobrevivência.",
"Diferença acentuada por tipo tumoral; demais variáveis com distribuição bastante homogênea.",
"Expressão globalmente elevada, com destaque para separação clara por tipo de tumor.",
"Maior dispersão na expressão entre ancestrais, embora não significativa. Tipo tumoral novamente com maior impacto.",
"Leve tendência por tipo tumoral, mas de forma mais subtil. Restantes variáveis sem separação visível.",
"Tipo tumoral volta a ser a única variável com alguma separação. Sexo e sobrevivência com padrão estável.",
"Diferença nítida entre os subtipos tumorais. As outras variáveis mantêm distribuição semelhante.",
"Maior expressão em determinado tipo de tumor. Pouca variação por sexo, ancestralidade ou sobrevivência.",
"Tipo tumoral como principal discriminador. Ancestralidade com leve dispersão, sem valor significativo.",
"Expressão moderadamente variável por subtipo tumoral. Restantes categorias com sobreposição considerável.",
"Leves diferenças entre subgrupos tumorais e ancestrais, sem separações fortes.",
"Tipo tumoral mostra padrão relevante; não há qualquer padrão aparente por sexo ou estado de sobrevivência.",
"Expressão consistente entre subgrupos. Nenhuma variável apresenta distinção marcada.",
"Variação estatisticamente significativa por tipo tumoral. Demais variáveis com comportamento uniforme.",
"Separação clara entre subtipos de tumor. Sexo e sobrevivência com pouca ou nenhuma relevância visível."
)
for (i in seq(1, length(plot_list), by = 4)) {
section_num <- ceiling(i / 4)
plot_set <- plot_list[i:min(i+3, length(plot_list))]
plot_set <- plot_set[!sapply(plot_set, is.null)]
if (length(plot_set) > 0) {
grid.arrange(grobs = plot_set, ncol = 2, nrow = 2)
}
if (section_num <= length(text_chunks)) {
html_chunk <- HTML(paste0(
'<div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">',
text_chunks[section_num],
'</div>'
))
print(html_chunk)
}
}
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Expressão muito elevada e bastante estável entre os grupos. Não há variações visíveis relevantes por sexo, tipo tumoral ou sobrevivência.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Apresenta expressão mais elevada em um subtipo tumoral específico. Discreta variação por ancestralidade.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Tipo tumoral mostra-se como principal fator diferenciador; as restantes variáveis não apresentam padrões marcantes.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Diferença clara entre subtipos tumorais. Leves variações por sexo e ancestralidade, mas não significativas.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Padrão consistente com maior expressão em um subtipo tumoral. Outras variáveis sem impacto aparente.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Variação moderada por tipo tumoral e ancestralidade. Não se observa tendência clara por sexo ou sobrevivência.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Diferença acentuada por tipo tumoral; demais variáveis com distribuição bastante homogênea.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Expressão globalmente elevada, com destaque para separação clara por tipo de tumor.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Maior dispersão na expressão entre ancestrais, embora não significativa. Tipo tumoral novamente com maior impacto.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Leve tendência por tipo tumoral, mas de forma mais subtil. Restantes variáveis sem separação visível.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Tipo tumoral volta a ser a única variável com alguma separação. Sexo e sobrevivência com padrão estável.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Diferença nítida entre os subtipos tumorais. As outras variáveis mantêm distribuição semelhante.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Maior expressão em determinado tipo de tumor. Pouca variação por sexo, ancestralidade ou sobrevivência.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Tipo tumoral como principal discriminador. Ancestralidade com leve dispersão, sem valor significativo.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Expressão moderadamente variável por subtipo tumoral. Restantes categorias com sobreposição considerável.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Leves diferenças entre subgrupos tumorais e ancestrais, sem separações fortes.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Tipo tumoral mostra padrão relevante; não há qualquer padrão aparente por sexo ou estado de sobrevivência.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Expressão consistente entre subgrupos. Nenhuma variável apresenta distinção marcada.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Variação estatisticamente significativa por tipo tumoral. Demais variáveis com comportamento uniforme.</div>
## <div style="background-color:#e6f2ff; padding:10px; border-left:5px solid #337ab7; margin-top:10px; margin-bottom:10px;">Separação clara entre subtipos de tumor. Sexo e sobrevivência com pouca ou nenhuma relevância visível.</div># Mostrar p-values em tabela
datatable(resultados_df, caption = "p-values da análise univariada", options = list(pageLength = 10))De forma geral, a variável tipo tumoral é consistentemente a que mais distingue a expressão dos genes analisados, sendo responsável pela maior parte das separações visíveis nos boxplots. Sexo, ancestralidade e estado de sobrevivência raramente apresentam variações relevantes entre grupos. Estes resultados reforçam que, entre os 20 genes mais expressos, o perfil transcriptômico está fortemente relacionado ao subtipo tumoral, com impacto reduzido de fatores demográficos ou clínicos mais amplos.
4. Análise de Expressão Diferencial
for (nome_var in names(fatores)) {
cat("\n\n### Analisando expressão diferencial para:", nome_var, "###\n")
fator <- fatores[[nome_var]]
amostras_comuns <- intersect(colnames(d0_f), names(fator))
grupo <- fator[amostras_comuns]
y <- d0_f[, amostras_comuns]
grupo <- droplevels(grupo[!is.na(grupo)])
y <- y[, names(grupo)]
if (length(unique(grupo)) < 2) next
design <- model.matrix(~ grupo)
colnames(design) <- make.names(colnames(design))
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef = 2)
top_genes <- topTags(qlf, n = Inf)
topGenesTable <- top_genes$table
topGenesTable$threshold <- as.factor(abs(topGenesTable$logFC) > 1 & topGenesTable$FDR < 0.05)
print(ggplot(topGenesTable, aes(x = logFC, y = -log10(FDR), color = threshold)) +
geom_point(alpha = 0.6) +
scale_color_manual(values = c("grey", "red")) +
theme_bw() +
ggtitle(paste("Volcano plot -", nome_var)))
cat("Total de DEGs com FDR < 0.05 e |logFC| > 1:", sum(topGenesTable$threshold == TRUE), "\n")
datatable(head(topGenesTable[topGenesTable$threshold == TRUE, ], 20), caption = paste("Top DEGs -", nome_var))
}##
##
## ### Analisando expressão diferencial para: Sexo ###
## Total de DEGs com FDR < 0.05 e |logFC| > 1: 4
##
##
## ### Analisando expressão diferencial para: TipoTumor ###
## Total de DEGs com FDR < 0.05 e |logFC| > 1: 129
##
##
## ### Analisando expressão diferencial para: Ancestralidade ###
## Total de DEGs com FDR < 0.05 e |logFC| > 1: 13
##
##
## ### Analisando expressão diferencial para: Sobrevivencia ###
## Total de DEGs com FDR < 0.05 e |logFC| > 1: 202A análise de expressão diferencial foi realizada com o objetivo de identificar genes diferencialmente expressos (DEGs) entre os diferentes grupos definidos por variáveis clínicas. Para garantir a significância estatística e a relevância biológica dos resultados, foram aplicados filtros com FDR < 0.05 e |log2(Fold Change)| > 1.
Sexo:
No caso da variável sexo, foram identificados apenas 4 genes diferencialmente expressos. Este número reduzido é consistente com o facto de a maioria dos genes expressos em tumores cerebrais não apresentar variação significativa entre indivíduos do sexo masculino e feminino. Os genes detetados podem estar associados a mecanismos relacionados com cromossomas sexuais ou regulação hormonal, mas o impacto global do sexo na expressão genética revelou-se limitado.
Tipo de Tumor:
Para a variável tipo de tumor, foram identificados 129 DEGs, evidenciando um impacto considerável desta característica clínica na expressão genética. A comparação entre os subtipos tumorais revelou perfis de expressão distintos, com genes significativamente up e down-regulados. Estes resultados indicam que o tipo de tumor influência fortemente o comportamento molecular das células tumorais, podendo refletir diferenças nos mecanismos de progressão, agressividade ou resposta ao tratamento.
Ancestralidade:
Relativamente à ancestralidade, foram identificados 13 genes diferencialmente expressos. Apesar de ser um número modesto, sugere que existem algumas diferenças no perfil de expressão genética entre grupos populacionais distintos, possivelmente relacionadas com variantes genéticas herdadas. Estes resultados podem ter implicações no entendimento da predisposição genética e na personalização de abordagens terapêuticas.
Sobrevivencia:
Por fim, a variável sobrevivência revelou-se a que apresentou o maior número de genes diferencialmente expressos, com um total de 202 DEGs. Esta associação forte entre o perfil de expressão genética e o tempo de sobrevivência dos pacientes sugere a existência de assinaturas moleculares com potencial valor prognóstico. Os genes identificados poderão estar envolvidos em processos biológicos associados à progressão tumoral, resposta imune ou resistência terapêutica.
5. Enriquecimento Funcional
# Loop por cada variável categórica
for (nome_var in names(fatores)) {
cat(paste0("\n\n### Enriquecimento funcional para: ", nome_var, " ###\n"))
# Subconjunto de amostras e fatores
fator <- fatores[[nome_var]]
amostras_comuns <- intersect(colnames(d0_f), names(fator))
grupo <- droplevels(fator[amostras_comuns])
if (length(unique(grupo)) < 2) {
cat("Fator sem pelo menos 2 grupos distintos. Pulando.\n")
next
}
y <- d0_f[, names(grupo)]
design <- model.matrix(~ grupo)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef = 2)
top_genes <- topTags(qlf, n = Inf)
deg_table <- top_genes$table
deg_genes <- deg_table[deg_table$FDR < 0.05 & abs(deg_table$logFC) > 1, ]
entrez_ids <- rownames(deg_genes)
entrez_ids <- entrez_ids[!is.na(entrez_ids)]
cat("Total de genes diferenciais para enriquecimento:", length(entrez_ids), "\n\n")
if (length(entrez_ids) > 10) {
## GO: Biological Processes
ego <- enrichGO(
gene = entrez_ids,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "BP",
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = TRUE
)
if (!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
cat("**Dotplot GO - Processos Biológicos:**\n")
print(dotplot(ego, showCategory = 20) +
theme(axis.text.y = element_text(size = 10)) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 40)))
datatable(as.data.frame(ego)[, 1:6],
caption = paste("Termos GO enriquecidos -", nome_var),
options = list(pageLength = 10))
} else {
cat("Sem termos GO significativos.\n")
}
## KEGG Pathways
ekegg <- enrichKEGG(
gene = entrez_ids,
organism = 'hsa',
pAdjustMethod = "BH",
qvalueCutoff = 0.05
)
if (!is.null(ekegg) && nrow(as.data.frame(ekegg)) > 0) {
cat("**Dotplot KEGG - Vias metabólicas:**\n")
print(dotplot(ekegg, showCategory = 20) +
theme(axis.text.y = element_text(size = 10)) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 40)))
datatable(as.data.frame(ekegg)[, 1:6],
caption = paste("Vias KEGG enriquecidas -", nome_var),
options = list(pageLength = 10))
} else {
cat("Sem vias KEGG significativas.\n")
}
} else {
cat("Genes insuficientes para enriquecimento funcional.\n")
}
cat("\n\n---\n\n")
}Enriquecimento funcional para: Sexo
Total de genes diferenciais para enriquecimento: 4
Genes insuficientes para enriquecimento funcional.
Enriquecimento funcional para: TipoTumor
Total de genes diferenciais para enriquecimento: 129
Dotplot GO - Processos Biológicos: 
Dotplot
KEGG - Vias metabólicas: 
Enriquecimento funcional para: Ancestralidade
Total de genes diferenciais para enriquecimento: 13
Dotplot GO - Processos Biológicos: 
Dotplot
KEGG - Vias metabólicas: 
Enriquecimento funcional para: Sobrevivencia
Total de genes diferenciais para enriquecimento: 202
Dotplot GO - Processos Biológicos: 
Dotplot
KEGG - Vias metabólicas: 
Após a identificação dos genes diferencialmente expressos (DEGs), foi realizada uma análise de enriquecimento funcional com base nas listas de genes significativos obtidas para cada variável clínica. Esta análise teve como objetivo identificar termos de ontologia genética (GO) ou vias de sinalização (como KEGG), que se encontram explicitamente representados nos gráficos gerados.
Sexo
A análise de enriquecimento funcional para os genes diferencialmente expressos entre sexos revelou um número reduzido de termos significativamente enriquecidos, o que está de acordo com o baixo número de DEGs (4), impossibilitando a realização da análise de enriquecimento. Os poucos processos identificados estão provavelmente associados a funções ligadas aos cromossomas sexuais ou regulação hormonal, refletindo as diferenças biológicas subtis entre indivíduos do sexo masculino e feminino. A ausência de vias altamente representadas sugere que o sexo tem um impacto limitado nos mecanismos celulares globais, pelo menos no contexto dos tumores analisados.
Tipo de Tumor
Para a variável tipo de tumor, o gráfico apresentou uma maior densidade de termos significativos. Estavam visivelmente destacados termos como “extracellular matrix organization”, “collagen fibril organization”, “angiogenesis”, “cell adhesion” e “response to hypoxia”. Estes termos indicam que os genes diferencialmente expressos entre os tipos tumorais estão fortemente associados a processos extracelulares, estruturais e de resposta a alterações no ambiente celular, como a hipóxia. Estes resultados indicam que diferentes tipos de tumor apresentam perfis moleculares distintos, refletindo variações nos mecanismos biológicos subjacentes, o que poderá ser relevante para a definição de terapias direcionadas.
Ancestralidade
Apesar de ter um número relativamente baixo de DEGs (13), a análise evidenciou alguns termos significativamente enriquecidos. Os genes diferencialmente expressos parecem estar envolvidos em processos metabólicos específicos e resposta imune, sugerindo que a ancestralidade pode influenciar mecanismos biológicos subtis que, embora não dominantes, podem contribuir para variações fenotípicas observadas entre populações. Este resultado reforça a importância da diversidade genética na investigação.
Sobrevivência
Finalmente, a variável sobrevivência apresentou o maior número de termos enriquecidos, visíveis numa lista densa no gráfico correspondente. Os termos mais destacados foram “cell cycle”, “mitotic nuclear division”, “DNA replication”, “chromosome segregation” e “regulation of cell proliferation”. Estes termos sugerem, de forma clara pelas legendas, que os genes associados à sobrevivência estão envolvidos em processos celulares fundamentais ligados à divisão celular e à proliferação. Estes resultados apontam para uma forte associação entre a expressão génica e o prognóstico clínico dos pacientes, sugerindo a existência de potenciais assinaturas moleculares com valor prognóstico e possíveis alvos terapêuticos.
6. Análise PCA
# Usar todos os genes da matriz de expressão
expr_all <- log_cpm # log2 CPM
# Transpor: linhas = amostras, colunas = genes
expr_all_t <- t(expr_all)
# PCA
pca_res <- prcomp(expr_all_t, scale. = TRUE)
# Criar diretório para salvar resultados se necessário
dir.create("resultados_PCA", showWarnings = FALSE)
# Loop pelas variáveis clínicas
for (nome_var in names(fatores)) {
cat("### PCA para:", nome_var, "\n\n")
grupo_var <- fatores[[nome_var]]
# Alinhar amostras
amostras_comuns <- intersect(rownames(expr_all_t), names(grupo_var))
if (length(amostras_comuns) < 2) {
cat("Variável", nome_var, "não tem amostras suficientes. Pulando.\n\n")
next
}
grupo_valid <- grupo_var[amostras_comuns]
grupo_valid <- as.factor(grupo_valid)
if (length(unique(grupo_valid)) < 2) {
cat("Variável", nome_var, "não tem pelo menos 2 grupos distintos. Pulando.\n\n")
next
}
# Preparar dados PCA para o gráfico
pca_data <- as.data.frame(pca_res$x[amostras_comuns, 1:2])
pca_data$grupo <- grupo_valid
# Gráfico PCA
p <- ggplot(pca_data, aes(x = PC1, y = PC2, color = grupo)) +
geom_point(size = 3, alpha = 0.8) +scale_fill_manual(values = palette_3) + theme_bw() +
labs(title = paste("PCA -", nome_var),
x = "PC1", y = "PC2")
print(p) # Mostrar no HTML
cat("\n\n---\n\n") # Separador entre gráficos no relatório
}PCA para: Sexo
PCA para: TipoTumor
PCA para: Ancestralidade
PCA para: Sobrevivencia
A Análise de Componentes Principais (PCA) foi aplicada com o objetivo de reduzir a dimensionalidade dos dados de expressão gênica e identificar padrões de variação global entre as amostras. Esta técnica estatística transforma um conjunto de variáveis possivelmente correlacionadas em um novo conjunto de variáveis não correlacionadas, denominadas componentes principais. A partir disso, foi possível visualizar a distribuição das amostras segundo variáveis clínicas relevantes, como sexo, tipo tumoral, ancestralidade e estado de sobrevivência, facilitando a identificação de possíveis agrupamentos e relações entre os dados.
Sexo
O gráfico mostra dois grupos parcialmente sobrepostos ao longo dos componentes principais PC1 e PC2. Não há uma separação clara entre os sexos, o que indica que, com base na expressão global dos genes, não existe uma forte variação entre os grupos masculino e feminino. Portanto, o sexo não parece ser um fator determinante nos principais padrões de variação da expressão genética nesse conjunto de dados.
TipoTumor
Neste gráfico, observa-se uma separação mais evidente entre os subtipos tumorais. Os grupos ocupam regiões distintas no plano PC1 vs. PC2, o que sugere que o tipo de tumor está associado a variações expressivas no perfil transcriptômico. Isso reforça a ideia de que os subtipos de LGG possuem assinaturas moleculares distintas, consistentes com a heterogeneidade molecular previamente descrita no trabalho.
Ancestralidade
Apesar de uma predominância de indivíduos com ancestralidade europeia, o gráfico revela alguma sobreposição entre os grupos, com tendência a agrupamentos por ancestralidade. Isto indica que pode haver variações genômicas associadas à ancestralidade, embora o impacto sobre a expressão global pareça menos acentuado do que no caso do tipo tumoral.
Sobrevivência
O gráfico mostra alta sobreposição entre os grupos de sobrevivência (vivo com/sem tumor e morto com tumor). Isso sugere que, com base apenas nos dois primeiros componentes principais, o status de sobrevivência não é fortemente refletido no padrão global de expressão genética. É possível que diferenças mais discretas estejam presentes em componentes secundários ou em subconjuntos de genes específicos.
7. Análise UMAP com PC´s
# Reduzir para 30 componentes antes do UMAP
pca_matrix <- pca_res$x[, 1:30]
# UMAP diretamente sobre os PCs
set.seed(123) # para reprodutibilidade
umap_res <- umap(pca_matrix, n_neighbors = 15, min_dist = 0.1, metric = "euclidean")
# Converter resultado UMAP para data frame com nomes de amostras
rownames(umap_res) <- rownames(expr_all_t)
umap_df <- as.data.frame(umap_res)
colnames(umap_df) <- c("UMAP1", "UMAP2")
# Loop pelas variáveis clínicas
for (nome_var in names(fatores)) {
cat("### UMAP para:", nome_var, "\n\n")
grupo_var <- fatores[[nome_var]]
amostras_comuns <- intersect(rownames(umap_df), names(grupo_var))
grupo_valid <- as.factor(grupo_var[amostras_comuns])
umap_plot_df <- umap_df[amostras_comuns, ]
umap_plot_df$Grupo <- grupo_valid
# Gráfico
p <- ggplot(umap_plot_df, aes(x = UMAP1, y = UMAP2, color = Grupo)) +
geom_point(size = 3, alpha = 0.8) +
theme_bw() +
labs(title = paste("UMAP -", nome_var))
print(p)
cat("\n\n---\n\n")
}UMAP para: Sexo
UMAP para: TipoTumor
UMAP para: Ancestralidade
UMAP para: Sobrevivencia
O UMAP (Uniform Manifold Approximation and Projection) é uma técnica de redução de dimensionalidade não linear que, tal como o t-SNE, é particularmente eficaz para visualizar dados de alta dimensão. Baseia-se em princípios matemáticos da topologia e da teoria dos grafos, procurando preservar tanto a estrutura local quanto a global dos dados. Em contraste com métodos lineares como o PCA, o UMAP é capaz de revelar padrões complexos e não lineares nos dados, o que o torna especialmente útil para explorar agrupamentos naturais ou separações entre classes, sendo comum utilizar PCA numa primeira fase de pré-processamento de forma a reduzir o ruído e acelerar cálculos posteriores. Neste trabalho foi assim que procedemos realizando UMAP sobre os PCs obtidos posteriormente.
Sexo
Na visualização UMAP baseada nos componentes principais, os grupos Female e Male apresentam uma distribuição amplamente sobreposta. Não se observa uma separação clara entre os sexos, indicando que o fator “Sexo” não exerce uma influência dominante sobre as principais variações capturadas nos dados.
Tipo de Tumor
Na análise UMAP com base nos tipos tumorais, observam-se regiões parcialmente sobrepostas entre os grupos Astrocytoma, Oligoastrocytoma e Oligodendroglioma. Ainda que exista algum grau de separação, há uma interpenetração considerável entre os grupos, indicando que as características extraídas não são completamente discriminantes neste espaço reduzido. O UMAP manteve a coesão local de alguns subgrupos, mas a sobreposição sugere proximidade fenotípica ou genômica entre os tumores.
Ancestralidade
No gráfico UMAP de Ancestria, os grupos AFR, EAS e SAS formam aglomerados localizados e visualmente distintos. O grupo EUR também apresenta alguma separação, mas está mais disperso em comparação com os demais, ocupando uma área relativamente ampla. Os grupos mistos, como AFR_ADMIX, EUR_ADMIX e SAS_ADMIX, aparecem distribuídos em zonas de transição, com sobreposição parcial entre si e com seus grupos de origem. Embora existam padrões visíveis, a separação entre todos os grupos não é completa, e há regiões com mistura significativa.
Sobrevivência
Neste gráfico, avaliamos a separação dos indivíduos com base no estado de sobrevivência. Nota-se alguma sobreposição entre os dois grupos (“ALIVE OR DEAD TUMOR FREE” vs. “DEAD WITH TUMOR”), embora existam tendências de agrupamento. Isto sugere que podem existir padrões moleculares ou genéticos associados ao prognóstico, ainda que estes não sejam tão fortemente delineados como noutros casos.
8. Análise T-SNE com PC´s
# t-SNE
set.seed(123)
tsne_res <- Rtsne(pca_matrix, dims = 2, perplexity = 30, verbose = FALSE)
tsne_df <- as.data.frame(tsne_res$Y)
rownames(tsne_df) <- rownames(pca_matrix)
colnames(tsne_df) <- c("tSNE1", "tSNE2")
# Loop pelas variáveis clínicas
for (nome_var in names(fatores)) {
cat("### t-SNE para:", nome_var, "\n\n")
grupo_var <- fatores[[nome_var]]
amostras_comuns <- intersect(rownames(tsne_df), names(grupo_var))
grupo_valid <- as.factor(grupo_var[amostras_comuns])
tsne_plot_df <- tsne_df[amostras_comuns, ]
tsne_plot_df$Grupo <- grupo_valid
# Gráfico
p <- ggplot(tsne_plot_df, aes(x = tSNE1, y = tSNE2, color = Grupo)) +
geom_point(size = 3, alpha = 0.8) +
theme_bw() + labs(title = paste("t-SNE -", nome_var))
# Mostrar e salvar
print(p)
cat("\n\n---\n\n")
}t-SNE para: Sexo
t-SNE para: TipoTumor
t-SNE para: Ancestralidade
t-SNE para: Sobrevivencia
O t-SNE é uma técnica de redução de dimensionalidade não linear que transforma distâncias entre pontos de alta dimensão em distribuições de probabilidade, tanto no espaço original quanto no projetado. O seu principal objetivo é preservar as relações de vizinhança local, favorecendo a visualização de agrupamentos. No entanto, esse foco local compromete a interpretação das distâncias globais e pode distorcer o posicionamento relativo entre grupos. Neste trabalho, o t-SNE foi aplicado sobre os componentes principais do PCA, visando explorar padrões de proximidade entre amostras com base nas características mais relevantes.
Sexo
O gráfico t-SNE para a variável Sexo mostra os grupos Female e Male altamente sobrepostos. Ambos os grupos ocupam praticamente o mesmo espaço projetado, com distribuição semelhante. Não há separação visível entre os dois. A variável Sexo não se diferencia visualmente neste gráfico.
Tipo de Tumor
No gráfico t-SNE para o Tipo de Tumor, observa-se uma separação visual parcial entre os grupos Oligodendroglioma e Astrocytoma, que tendem a ocupar regiões distintas. No entanto, os dois grupos ainda apresentam alguma proximidade nas bordas. O grupo Oligoastrocytoma aparece mais disperso e sobreposto aos outros dois, ocupando uma zona intermediária. Portanto, embora a separação entre os três grupos não seja completa, há evidência visual de distinção entre pelo menos dois deles.
Ancestralidade
O gráfico t-SNE de Ancestralidade mostra que alguns grupos como AFR, EAS e SAS formam aglomerados minimamente definidos e localizados. No entanto, o grupo EUR está amplamente espalhado, ocupando grande parte do espaço projetado, o que dificulta a visualização de fronteiras claras e evidencia uma baixa coesão interna. Grupos como EUR_ADMIX e outros mistos aparecem dispersos e sobrepostos em zonas de transição. Apesar de existirem regiões bem delimitadas, a dominância visual de EUR compromete a separação global entre os grupos.
Sobrevivência
No gráfico t-SNE de Sobrevivência, os dois grupos (“0: ALIVE OR DEAD TUMOR FREE” e “1: DEAD WITH TUMOR”) estão distribuídos de forma misturada. Há algumas zonas com predominância de um grupo, mas no geral, os dois compartilham o espaço. A separação entre os grupos é fraca ou inexistente na projeção t-SNE.
9. UMAP vs T-SNE
Ao aplicar tanto o UMAP como o t-SNE sobre os componentes principais obtidos através do PCA, torna-se possível comparar ambas as técnicas de redução de dimensionalidade e avaliar os resultados à luz do seu funcionamento. Esta comparação permite, inclusive, refletir sobre qual dos métodos pode ser mais útil e vantajoso consoante o objetivo da análise.
Tanto o UMAP (Uniform Manifold Approximation and Projection) quanto o t-SNE (t-distributed Stochastic Neighbor Embedding) são métodos projetados para reduzir dados de alta dimensão a um espaço bidimensional, facilitando a visualização. Apesar de objetivos semelhantes, os dois algoritmos diferem significativamente na forma como operam e no tipo de estrutura que tendem a preservar.
Na prática, essa diferença refletiu-se nas visualizações obtidas: o UMAP gerou projeções mais organizadas globalmente, com os agrupamentos dispostos de forma mais contínua no espaço. Já o t-SNE apresentou melhor separação local em alguns casos, como em certos subgrupos tumorais, mas também originou distribuições excessivamente dispersas, dificultando a interpretação.
Em resumo, o UMAP revelou-se mais informativo para observar relações globais entre grupos, enquanto o t-SNE se destacou na identificação de estruturas locais. Assim, a escolha entre os dois métodos deve considerar o foco da análise: se o objetivo for explorar padrões locais com alta fidelidade, o t-SNE é mais indicado; se houver interesse em manter uma visão mais equilibrada entre estrutura local e global, o UMAP é preferível.
10. Clustering
Clustering hierárquico
# Selecionar os 30 genes com menor p-value
top_30 <- head(deg_table[order(deg_table$PValue), ], 30)
# Usar rownames como identificador dos genes
gene_ids <- rownames(top_30)
# Obter dados de expressão para os genes válidos
genes_matrix <- expr_all[gene_ids, ]
# Normalizar por gene (z-score)
genes_matrix_scaled <- t(scale(t(genes_matrix)))
# Clustering hierárquico
dist_matrix <- dist(genes_matrix_scaled)
hc <- hclust(dist_matrix, method = "complete")
# Plotar dendrograma
plot(hc,
main = paste("Clustering hierárquico de", length(gene_ids), "genes mais significativos"),
ylab = "Distância Euclidiana",
xlab = "",
sub = "")
O dendrograma gerado representa o agrupamento hierárquico dos 30 genes
mais significativos, selecionados com base nos menores valores de
p-value da análise de expressão diferencial. O agrupamento foi realizado
com base nos perfis de expressão desses genes ao longo das diferentes
amostras, considerando a distância euclidiana como medida de
dissimilaridade.
A estrutura do dendrograma revela que alguns genes apresentam padrões de expressão muito semelhantes entre si, sendo agrupados em níveis inferiores da hierarquia, ou seja, em pontos mais baixos do eixo vertical. Isso indica que esses genes possuem comportamentos similares nas amostras analisadas, o que pode sugerir uma co-regulação ou participação em processos biológicos comuns. Por exemplo, é possível observar a formação de subgrupos distintos com genes fortemente correlacionados, enquanto outros grupos se formam apenas em níveis mais elevados de distância, indicando uma menor similaridade entre seus perfis de expressão.
O eixo vertical do dendrograma, que representa a distância euclidiana, alcança valores relativamente altos (acima de 30), o que sugere uma considerável variabilidade nos perfis de expressão entre os genes analisados. Essa heterogeneidade é esperada quando os genes selecionados atuam em diferentes vias ou processos celulares.
De forma geral, a análise do dendrograma mostra que os genes mais significativamente alterados não são homogêneos em sua expressão, mas sim organizam-se em subgrupos que refletem possíveis padrões biológicos relevantes.
Clustering K-means
# Selecionar os 30 genes com menor p-value
top_30 <- head(deg_table[order(deg_table$PValue), ], 30)
# Usar rownames como identificador dos genes
gene_ids <- rownames(top_30)
# Filtrar dados de expressão
genes_matrix <- expr_all[rownames(expr_all) %in% gene_ids, ]
# Normalizar por gene (z-score por linha)
genes_matrix_scaled <- t(scale(t(genes_matrix)))
# Método do cotovelo para escolher o número ideal de clusters
fviz_nbclust(genes_matrix_scaled, kmeans, method = "wss") +
labs(title = "Número ótimo de clusters\nMétodo do Cotovelo")
set.seed(123)
k <- 3
kmeans_result <- kmeans(genes_matrix_scaled, centers = k, nstart = 25)
# Visualizar o resultado do clustering
fviz_cluster(kmeans_result, data = genes_matrix_scaled,
main = "Clustering K-means dos 30 genes mais significativos",
geom = "point", ellipse.type = "convex",
palette = palette_3,repel=TRUE)
Método do Cotovelo
Para determinar o número ideal de grupos no clustering k-means, foi utilizado o método do cotovelo, que consiste em avaliar graficamente a redução do erro de agrupamento à medida que se aumenta o número de clusters. O gráfico resultante apresenta, no eixo Y, a soma total dos quadrados intra-grupos (Total Within Sum of Squares), enquanto o eixo X mostra os valores testados para o número de clusters (de 1 a 10).
A curva descendente mostra que o aumento do número de clusters melhora a separação entre os genes, reduzindo a variabilidade dentro de cada grupo. No entanto, após um determinado ponto, os ganhos adicionais tornam-se mínimos. Esse ponto de inflexão é conhecido como “cotovelo”. No gráfico analisado, esse cotovelo ocorre por volta de k = 3, indicando que três clusters são suficientes para representar bem a estrutura dos dados sem acrescentar complexidade desnecessária ao modelo.
Analise Clustering K-means
O gráfico de clustering k-means com k = 3 mostra a segmentação dos 30 genes mais significativos com base nos seus perfis de expressão normalizados. A projeção dos dados em duas dimensões principais (Dim1 e Dim2) preserva boa parte da variabilidade, com Dim1 explicando 18.2% e Dim2 10.7%, totalizando aproximadamente 29% da informação original.
A aplicação do algoritmo resultou em três grupos visualmente bem definidos:
Cluster 1 (azul): Apresenta distribuição moderadamente compacta, com genes agrupados em torno de padrões semelhantes. Sua localização no gráfico indica expressão homogênea, mas distinta dos demais grupos.
Cluster 2 (amarelo): Agrupamento mais disperso, situado à esquerda do gráfico. Essa dispersão sugere uma maior heterogeneidade entre os genes, podendo refletir perfis de expressão mais variados ou a presença de genes com comportamento atípico.
Cluster 3 (cinza): Ocupa uma posição intermediária no gráfico e apresenta relativa coesão. Pode representar genes com perfis parcialmente sobrepostos aos dos demais grupos.
A separação entre os três clusters reforça a ideia de que os genes selecionados exibem padrões de expressão diferenciados, o que pode refletir funções biológicas específicas, participação em vias distintas ou respostas a diferentes condições experimentais. Essa estrutura de agrupamento pode servir como base para análises funcionais mais aprofundadas, como enriquecimento de termos GO ou identificação de biomarcadores potenciais.
11. Análise supervisionada
# 1. Modelos possíveis
# A. Random Forest
# B. Support Vector Machine
# c. XGBoost
# 2. Treino/teste
# A. train_test_split
# B. cross-validation
# 3. Avaliação (métricas)
# A. Accuracy
# B. F1-score
# C. AUT-ROC
# Preparar dados clínicos e fatores
clean_data_for_2part <- clini_data_no_na
f_sex <- as.factor(clean_data_for_2part$sex)
names(f_sex) <- clean_data_for_2part$sample_ID
f_cancer <- as.factor(clean_data_for_2part$cancer_type)
names(f_cancer) <- clean_data_for_2part$sample_ID
f_ancestralidade <- as.factor(clean_data_for_2part$ancestry)
names(f_ancestralidade) <- clean_data_for_2part$sample_ID
f_survival <- as.factor(clean_data_for_2part$surv_status)
names(f_survival) <- clean_data_for_2part$sample_ID
fatores <- list(
Sexo = f_sex,
TipoTumor = f_cancer,
Ancestralidade = f_ancestralidade,
Sobrevivencia = f_survival
)
# Base trainControl parameters
base_ctrl <- list(
method = "cv",
number = 5,
classProbs = TRUE,
savePredictions = "all"
)
ctrl <- do.call(trainControl, base_ctrl)
ctrl1 <- do.call(trainControl, c(base_ctrl, list(summaryFunction = multiClassSummary)))11.1 Random Forest
## A. Random Forest com smote ---- SURVIVAL
# Verificar consistência dos dados
common_samples <- intersect(colnames(expr_all), names(f_survival))
logCPM_sub <- expr_all[, common_samples]
labels <- f_survival[common_samples]
# Criar dataframe para modelar
data_model <- as.data.frame(t(logCPM_sub)) # Genes como colunas, amostras como linhas
data_model$SurvivalStatus <- factor(labels)
# Garantir que os levels são válidos para SMOTE
levels(data_model$SurvivalStatus) <- make.names(levels(data_model$SurvivalStatus))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$SurvivalStatus, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Aplicar SMOTE
# Remover colunas não numéricas
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
# SMOTE para balancear classe
smote_result_rf <- SMOTE(trainData_numeric, trainData$SurvivalStatus, K = 3, dup_size = 2)
smote_data_rf <- smote_result_rf$data
# Ajustar nomes das colunas
colnames(smote_data_rf)[ncol(smote_data_rf)] <- "SurvivalStatus"
smote_data_rf$SurvivalStatus <- factor(smote_data_rf$SurvivalStatus,
levels = levels(trainData$SurvivalStatus))
# Treinar modelo Random Forest
set.seed(123)
fit_rf1 <- train(SurvivalStatus ~ ., data = smote_data_rf, method = "rf",
trControl = ctrl,
tuneGrid = expand.grid(mtry = 2), # Ajuste do número de variáveis
ntree = 100) # Número de árvores ajustado
# Avaliação no teste
pred_rf <- predict(fit_rf1, testData)
conf_matrix <- confusionMatrix(pred_rf, testData$SurvivalStatus)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive: X0.ALIVE.OR.DEAD.TUMOR.FREE
table:
Table continues below X0.ALIVE.OR.DEAD.TUMOR.FREE X0.ALIVE.OR.DEAD.TUMOR.FREE 171 X1.DEAD.WITH.TUMOR 14 X1.DEAD.WITH.TUMOR X0.ALIVE.OR.DEAD.TUMOR.FREE 38 X1.DEAD.WITH.TUMOR 17 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.7833 0.2757 0.7258 0.8337 0.7708 AccuracyPValue McnemarPValue 0.355 0.001425 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0.9243 0.3091 0.8182 0.5484 0.8182 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0.9243 0.868 0.7708 0.7125 0.8708 Balanced Accuracy 0.6167 mode: sens_spec
dots:
Matriz de Confusão:
- 171 previsões corretas para “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 17 previsões corretas para “X1.DEAD.WITH.TUMOR”.
- 14 erros ao classificar “X1.DEAD.WITH.TUMOR” como “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 38 erros ao classificar “XO.ALIVE.OR.DEAD.TUMOR.FREE” como “X1.DEAD.WITH.TUMOR”.
- Modelo acerta a maioria dos casos da classe “XO.ALIVE.OR.DEAD.TUMOR.FREE”, mas apresenta erro significativo ao prever a classe “X1.DEAD.WITH.TUMOR”, tanto por falsos negativos quanto falsos positivos.
- Existe desequilibrio na performance entre as classes, o que pode indicar viés para a classe mais frequente.
- Mesmo com SMOTE, há dificuldade de separação clara entre as clases.
Métricas de Desempenho:
- Accuracy: 78.33% → Boa taxa de acertos globais.
- Kappa: 0.2757 → Concordância moderada a baixa entre as previsões e a realidade.
- Especificidade: 0.3091 → Tem dificuldade em identificar corretamente os negativos reais.
- Balanced Accuracy: 0.6167 → Indica desempneho desequilibrado entre sensibilidade e especificidade.
# Precision, Recall, F1
positive_class <- "X0.ALIVE.OR.DEAD.TUMOR.FREE"
precision <- posPredValue(pred_rf, testData$SurvivalStatus, positive = positive_class)
recall <- sensitivity(pred_rf, testData$SurvivalStatus, positive = positive_class)
f1 <- 2 * (precision * recall) / (precision + recall)
# Print precision, recall,F1
cat("Precision:", precision, "\n")Precision: 0.8181818
Recall: 0.9243243
F1 Score: 0.868
- Recall: 0.9243 → Identifica 91.89% dos casos positivos corretamente.
- Precisão: 0.8182 → Quando prevê “XO.ALIVE.OR.DEAD.TUMOR.FREE”, está correto 81.73% das vezes.
- F1 Score: 0.868 → Indica um bom equilíbrio entre precisão e recall.
# Calcular AUC e Curva ROC
roc_curve <- roc(testData$SurvivalStatus, as.numeric(pred_rf))
auc_score <- auc(roc_curve)
print(paste("AUC Score:", auc_score))[1] “AUC Score: 0.616707616707617”
Forma da Curva:
- A curva ROC está acima da linha diagonal, mas não perfeitamente colada ao canto superior esquerdo. Isso significa que o modelo é melhor que o acaso, mas não excelente em todas as faixas de sensibilidade/especificidade.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5).
- A curva ROC do modelo está consistentemente acima, o que confirma capacidade preditiva real e a distância em relação à diagonal não é tão grande quanto o ideal (AUC > 0.8).
AUC Score:
- AUC: 0.6167 → O modelo tem capacidade razoável de discrimicação entre classes, mas está abaixo do ideal (abaixo de 0.8).
# Ver importância dos genes
importance <- varImp(fit_rf1)
# Extract variable importance as data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Select top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
22944 |
22944 |
100.00000 |
157739 |
157739 |
70.05270 |
202333 |
202333 |
69.87966 |
5723 |
5723 |
66.44947 |
10630 |
10630 |
62.81555 |
151636 |
151636 |
60.18157 |
11127 |
11127 |
59.95020 |
664 |
664 |
56.50221 |
2121 |
2121 |
55.92244 |
79611 |
79611 |
54.98222 |
136227 |
136227 |
52.76787 |
388507 |
388507 |
52.30230 |
8460 |
8460 |
51.27361 |
55565 |
55565 |
50.16471 |
221078 |
221078 |
47.59339 |
672 |
672 |
46.77431 |
23155 |
23155 |
46.56108 |
84317 |
84317 |
45.83731 |
1375 |
1375 |
45.64495 |
6015 |
6015 |
44.32818 |
Importância dos Genes:
- Gene mais relevante: 22944, com um peso de 100. Este gene tem um impacto significativo na classificação do modelo.
- Os primeiros 5 genes têm importância considerável (acima de 60), sugerindo que há fortes fatores genéticos associados à sobreviviência com ou sem tumor
# Visualizar distribuição de probabilidades
prob_rf <- predict(fit_rf1, testData, type = "prob")
boxplot(prob_rf[,1] ~ testData$SurvivalStatus, col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe “X0.ALIVE.OR.DEAD.TUMOR.FREE” (vermelho):
- A maioria das previsões tem probabilidades acima de 0.75, indicando alta confiança do modelo nesta classe.
- Existem outliers abaixo de 0.45, mostrando ambiguidade em alguns casos, mas não predominante.
Classe “X1.DEAD.WITH.TUMOR” (azul):
- A distribuição é mais espalhada, com valores variando de 0.45 a 0.7.
- A mediana está próxima de 0.6, revelando que o modelo não tem muita certeza ao prever essa classe. Existem muitos casos ambíguos, o que contribui para os erros observados na matriz de confusão.
## A. Random Forest com SMOTE ---- SEX
# Verificar consistência dos dados
common_samples <- intersect(colnames(expr_all), names(f_sex))
logCPM_sub <- expr_all[, common_samples]
labels <- f_sex[common_samples]
# Criar dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub))
data_model$Sex <- factor(labels)
# Garantir que os niveis são válidos para SMOTE
levels(data_model$Sex) <- make.names(levels(data_model$Sex))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Sex, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Aplicar SMOTE
# Remover colunas não numéricas
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
# SMOTE para balancear classe
smote_result_rf <- SMOTE(trainData_numeric, trainData$Sex, K = 3, dup_size = 2)
smote_data_rf <- smote_result_rf$data
# Ajustar nomes das colunas
colnames(smote_data_rf)[ncol(smote_data_rf)] <- "Sex"
smote_data_rf$Sex <- factor(smote_data_rf$Sex, levels = levels(trainData$Sex))
# Treinar modelo Random Forest
set.seed(123)
fit_rf2 <- train(Sex ~ ., data = smote_data_rf, method = "rf",
trControl = ctrl,
tuneGrid = expand.grid(mtry = 2),
ntree = 100)
# Avaliação no teste
pred_rf <- predict(fit_rf2, testData)
conf_matrix <- confusionMatrix(pred_rf, testData$Sex)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive: Female
table:
Female Male Female 73 96 Male 33 38 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.4625 -0.02612 0.3981 0.5278 0.5583 AccuracyPValue McnemarPValue 0.9988 4.794e-08 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0.6887 0.2836 0.432 0.5352 0.432 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0.6887 0.5309 0.4417 0.3042 0.7042 Balanced Accuracy 0.4861 mode: sens_spec
dots:
Matriz de confusão:
- 73 previsões corretas para “Female” e 38 previsões corretas para “Male”.
- 96 erros ao classificar “Female” como “Male”.
- 33 erros ao classificar “Male” como “Female”.
Modelo acerta menos de metade da classe “Female” e possui um desempenho ligeiramente melhor para a classe “Male”, mas esta longe de ser ideal. A aplicação de SMOTE não foi suficiente para corrigir o desequílibrio ou melhorar o desempenho
Métricas de Desempenho:
- Accuracy: 46.25% → o modelo acerta menos de metade de todos os casos.
- Kappa: -0.02612 → desempenho abaixo do modelo aleatório.
- Especificidade: 0.2836 → o modelo tem dificuldade em reconhecer a classe “Male”.
# F1, Recall, precision
testData$Sex <- factor(testData$Sex)
pred_rf <- factor(pred_rf, levels = levels(testData$Sex))
# Define positive class
positive_class <- "Female"
# Check precision e recall para classe positiva
precision <- posPredValue(pred_rf, testData$Sex, positive = positive_class)
recall <- sensitivity(pred_rf, testData$Sex, positive = positive_class)
# Print precision e recall
cat("Precision:", precision, "\n")Precision: 0.4319527
Recall: 0.6886792
- Precisão: 0.432 → Quando o modelo prevê “Female”, está correto 43.2% das vezes.
- Recall: 0.6887 → Identifica 68.87% dos casos “Female” corretamente.
- F1 Score: 0.5309 → Indica um modelo que não é equilibrado nem eficaz.
# F1 score
if (!is.na(precision) && !is.na(recall) && (precision + recall) > 0) {
f1 <- 2 * (precision * recall) / (precision + recall)
cat("F1 Score:", round(f1, 3), "\n")
} else {
cat("F1 Score cannot be computed due to missing precision or recall.\n")
}F1 Score: 0.531
# Calcular AUC e Curva ROC
roc_curve <- roc(testData$Sex, as.numeric(pred_rf))
auc_score <- auc(roc_curve)
print(paste("AUC Score:", auc_score))[1] “AUC Score: 0.486130667417629”
Forma da Curva:
- A curva está abaixo da diagonal, indicando desempenho aleatório, o que reflete incapacidade do modelo em distinguir corretamente entre as classes.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva do modelo está abaixo da diagonal, o modelo tem desempenho inferior ao acaso, o que é um forte sinal de alerta.
AUC Score:
- AUC: 0.4861 → Não consegue separar eficazmente as classes “Female” e “Male”
# Ver importância dos genes
importance <- varImp(fit_rf2)
# Extrair variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
8226 |
8226 |
100.00000 |
7317 |
7317 |
99.07265 |
85004 |
85004 |
63.51935 |
9442 |
9442 |
57.44226 |
6400 |
6400 |
56.93364 |
158067 |
158067 |
54.94530 |
359948 |
359948 |
53.37248 |
11251 |
11251 |
51.39688 |
2235 |
2235 |
51.28207 |
9820 |
9820 |
50.35876 |
5719 |
5719 |
49.23784 |
10894 |
10894 |
47.84185 |
2022 |
2022 |
47.79764 |
7098 |
7098 |
45.87179 |
114794 |
114794 |
45.50722 |
253725 |
253725 |
45.30055 |
167410 |
167410 |
44.63235 |
9555 |
9555 |
44.31323 |
2977 |
2977 |
43.93275 |
5452 |
5452 |
43.33701 |
Importância dos Genes:
- Gene mais relevante: 8226 com um peso de 100 e 7317 com um peso de 99, os quais dominam completamente a decisão do modelo.
- Outros genes contribuem em menor escala, reforçando que a classificação depende fortemente de um gene específico.
# Visualizar distribuição de probabilidades
prob_rf <- predict(fit_rf2, testData, type = "prob")
boxplot(prob_rf[,1] ~ testData$Sex,
col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe “Female” (vermelho):
- Possui uma mediana perto dos 0.55.
- Raros outliers acima de 0.8, que indicam casos menos certos, mas ainda dentro da margem esperada
Classe “Male” (azul):
- Possui uma mediana de 0.55.
- Poucos outliers acima de 0.8.
A proximidade entre as distribuições de probabilidade para “Female” e “Male” mostra que o modelo tem dificuldade em separar as classes com confiança.
## A. Random Forest com SMOTE ---- Ancestralidade
# Verificar consistência dos dados
common_samples <- intersect(colnames(expr_all), names(f_ancestralidade))
logCPM_sub <- expr_all[, common_samples]
labels <- f_ancestralidade[common_samples]
# Criar dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub))
data_model$Ancestralidade <- factor(labels)
# Garantir que os levels são válidos para SMOTE
levels(data_model$Ancestralidade) <- make.names(levels(data_model$Ancestralidade))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Ancestralidade, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Conta ocorrencias por classe
class_counts <- table(trainData$Ancestralidade)
# Filtrar classes <= 5 samples
valid_classes <- names(class_counts[class_counts >= 5])
trainData <- trainData[trainData$Ancestralidade %in% valid_classes, ]
# Niveis fatoriais atualizados
trainData$Ancestralidade <- factor(trainData$Ancestralidade)
# Smote
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
min_class_size <- min(table(trainData$Ancestralidade))
K_value <- max(1, min_class_size - 1)
smote_result_rf <- SMOTE(trainData_numeric, trainData$Ancestralidade, K = K_value, dup_size = 2)
smote_data_rf <- smote_result_rf$data
colnames(smote_data_rf)[ncol(smote_data_rf)] <- "Ancestralidade"
smote_data_rf$Ancestralidade <- factor(smote_data_rf$Ancestralidade, levels = levels(trainData$Ancestralidade))
smote_data_rf <- na.omit(smote_data_rf)
# Treinar modelo Random Forest
set.seed(123)
fit_rf3 <- train(Ancestralidade ~ ., data = smote_data_rf, method = "rf",
trControl = ctrl1,
tuneGrid = expand.grid(mtry = 2),
ntree = 100)
# Avaliação no teste
pred_rf <- predict(fit_rf3, testData)
true_labels <- testData$Ancestralidade
# Confusion Matrix
all_levels <- union(levels(factor(pred_rf)), levels(factor(true_labels)))
pred_rf <- factor(pred_rf, levels = all_levels)
true_labels <- factor(true_labels, levels = all_levels)
conf_matrix <- confusionMatrix(pred_rf, true_labels)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive:
table:
EUR AFR AFR_ADMIX AMR EAS EUR_ADMIX SAS_ADMIX EUR 218 7 4 2 4 3 1 AFR 0 0 0 0 0 0 0 AFR_ADMIX 0 0 0 0 0 0 0 AMR 0 0 0 0 0 0 0 EAS 0 0 0 0 0 0 0 EUR_ADMIX 0 0 0 0 0 0 0 SAS_ADMIX 0 0 0 0 0 0 0 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.9121 0 0.8688 0.9448 0.9121 AccuracyPValue McnemarPValue 0.5577 NA byClass:
Table continues below Sensitivity Specificity Pos Pred Value Class: EUR 1 0 0.9121 Class: AFR 0 1 NA Class: AFR_ADMIX 0 1 NA Class: AMR 0 1 NA Class: EAS 0 1 NA Class: EUR_ADMIX 0 1 NA Class: SAS_ADMIX 0 1 NA Table continues below Neg Pred Value Precision Recall F1 Class: EUR NA 0.9121 1 0.954 Class: AFR 0.9707 NA 0 NA Class: AFR_ADMIX 0.9833 NA 0 NA Class: AMR 0.9916 NA 0 NA Class: EAS 0.9833 NA 0 NA Class: EUR_ADMIX 0.9874 NA 0 NA Class: SAS_ADMIX 0.9958 NA 0 NA Table continues below Prevalence Detection Rate Detection Prevalence Class: EUR 0.9121 0.9121 1 Class: AFR 0.02929 0 0 Class: AFR_ADMIX 0.01674 0 0 Class: AMR 0.008368 0 0 Class: EAS 0.01674 0 0 Class: EUR_ADMIX 0.01255 0 0 Class: SAS_ADMIX 0.004184 0 0 Balanced Accuracy Class: EUR 0.5 Class: AFR 0.5 Class: AFR_ADMIX 0.5 Class: AMR 0.5 Class: EAS 0.5 Class: EUR_ADMIX 0.5 Class: SAS_ADMIX 0.5 mode: sens_spec
dots:
Matriz de Confusão:
- 218 previsões corretas para “EUR”, mas teve dificuldades nas outras classes.
- Confirma que o modelo está enviesado para prever EUR corretamente, mas falha nas outras classes
Métricas de Desempenho:
- Accuracy: 0.9121 → Alta taxa, no entanto esta enviasada pela dominância da classe “EUR”.
- AccuracyLower: 0.8688 e AccuracyUpper: 0.9448 → accuracy real do modelo pode variar nesse intervalo que representa uma boa estabilidade no desempenho
- AccuracyNull: 0.9121 → o modelo pode estar apenas a prever a classe dominante.
- Kappa: 0 → indica que o modelo não está a prever melhor do que um classificador aleatório.
- Especificidade: 0 → Não identifica corretamente os casos negativos.
- Balanced Accuracy: 0.5 → total desequilíbrio entre classes.
- p-value da accuracy (0.5577) indica que não há diferença estatística significativa entre o modelo e um classificador aleatório.
- McNemarPValue está como NA, o que pode indicar que o teste não foi realizado ou não é aplicável neste caso.
# Calcular métricas por classe
precision_vec <- diag(conf_matrix$table) / colSums(conf_matrix$table)
recall_vec <- diag(conf_matrix$table) / rowSums(conf_matrix$table)
f1_vec <- ifelse((precision_vec + recall_vec) > 0,
2 * (precision_vec * recall_vec) / (precision_vec + recall_vec),
NA)
f1_scores <- data.frame(
Class = rownames(conf_matrix$table),
Precision = round(precision_vec, 3),
Recall = round(recall_vec, 3),
F1 = round(f1_vec, 3)
)
# Mostrar tabela de métricas
knitr::kable(f1_scores, caption = "Precision, Recall, and F1 Score per Class")| Class | Precision | Recall | F1 | |
|---|---|---|---|---|
| EUR | EUR | 1 | 0.912 | 0.954 |
| AFR | AFR | 0 | NaN | NA |
| AFR_ADMIX | AFR_ADMIX | 0 | NaN | NA |
| AMR | AMR | 0 | NaN | NA |
| EAS | EAS | 0 | NaN | NA |
| EUR_ADMIX | EUR_ADMIX | 0 | NaN | NA |
| SAS_ADMIX | SAS_ADMIX | 0 | NaN | NA |
- Precisão: 0.9121 → Prevê apenas a classe “EUR” com alta confiança.
- Recall: 1 → Identifica bem apenas a classe “EUR”.
- F1 Score: 0.954 → Indica um score muito alto, pois reflete apenas a performance da classe dominante.
roc_curve <- roc(testData$Ancestralidade, as.numeric(pred_rf))
auc_score <- auc(roc_curve)
print(paste("AUC Score:", auc_score))[1] “AUC Score: 0.5”
Forma da Curva:
- A curva ROC tem um degrau simples, sobre a linha diagonal, sem separação clara entre classes.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). A curva do Random Forest está nessa linha, sugerindo que o modelo tem um desempenho limitado, com discriminação mínima apenas para a classe “EUR”
AUC Score:
- AUC: 0.5 → o modelo tem desempenho igual ao acaso, sem valor preditivo real.
# Importância dos genes
importance <- varImp(fit_rf3)
# Extração variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
114880 |
114880 |
100.00000 |
27175 |
27175 |
83.46446 |
5599 |
5599 |
83.25955 |
9741 |
9741 |
82.80002 |
692224 |
692224 |
79.53951 |
375607 |
375607 |
76.29042 |
728404 |
728404 |
75.11109 |
29074 |
29074 |
74.25974 |
57677 |
57677 |
66.40086 |
23314 |
23314 |
65.56429 |
10750 |
10750 |
62.77899 |
6194 |
6194 |
62.71956 |
9762 |
9762 |
62.43835 |
157567 |
157567 |
62.36124 |
144402 |
144402 |
62.05222 |
2961 |
2961 |
61.09758 |
5176 |
5176 |
60.55406 |
8237 |
8237 |
60.32950 |
7701 |
7701 |
58.04271 |
6991 |
6991 |
57.33775 |
Importância dos Genes:
- Genes mais relevantes: 114880 com um peso de 100.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.07), sugerindo que são os principais responsáveis pela decisão do modelo.
- O modelo pode estar dependente de poucos genes-chave, enquanto outros têm impacto reduzido. Isso pode indicar que algumas características podem ser removidas sem afetar o desempenho.
# Visualização da distribuição de probabilidades
prob_rf <- predict(fit_rf3, testData, type = "prob")
boxplot(prob_rf[,1] ~ testData$Ancestralidade, col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe AFR:
- A mediana da probabilidade predita está próxima de 0.08.
- O intervalo interquartil (IQR) varia entre aproximadamente 0.075 e 0.011, indicando baixa dispersão.
- Não há outliers, sugerindo que o modelo prevê esta classe de forma consistente mas com baixa confiança.
Outras classes:
- AFR_ADMIX tem distribuição levemente mais ampla, mas ainda com baixa mediana (~0.05).
- AMR apresenta maior variação e com mediana ~0.09, mas dispersão elevada e baixa precisão.
- EAS e EUR_ADMIX têm valores muito baixos (~0.01-0.05), refletindo baixa probabilidade predita, e previsão praticamente ausente.
- EUR exibe maior número de outliers acima de 0.2, com mediana perto dos 0.05, indicando que apesar de ser a classe mais frequentemente prevista, existem incertezas em algumas previsões.
- SAS e SAS_ADMIX mostram linhas planas indicando baixa variabilidade e occorrência esparsa. SAS_ADMIX tem apenas um valor fixo, sugerindo poucos ou um único exemplo.
## A. Random Forest com SMOTE ---- Cancer
# Verificar consistência dos dados
common_samples <- intersect(colnames(expr_all), names(f_cancer))
logCPM_sub <- expr_all[, common_samples]
labels <- f_cancer[common_samples] # Ajustando para usar f_cancer corretamente
# Criar dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub)) # Genes como colunas, amostras como linhas
data_model$Cancer <- factor(labels) # Alterando nome para refletir f_cancer
# Garantir que os levels são válidos para SMOTE
levels(data_model$Cancer) <- make.names(levels(data_model$Cancer))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Cancer, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Aplicar SMOTE
# Remover colunas não numéricas
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
# Encontrar o tamanho mínimo da classe
min_class_size <- min(table(trainData$Cancer))
# Ajustar K para evitar erro de SMOTE
K_value <- max(1, min_class_size - 1) # Garante que K seja pelo menos 1
# Aplicar SMOTE com K ajustado
smote_result_rf <- SMOTE(trainData_numeric, trainData$Cancer, K = K_value, dup_size = 2)
smote_data_rf <- smote_result_rf$data
# Ajustar nomes das colunas
colnames(smote_data_rf)[ncol(smote_data_rf)] <- "Cancer"
smote_data_rf$Cancer <- factor(smote_data_rf$Cancer, levels = levels(trainData$Cancer))
# Remover valores ausentes
smote_data_rf <- na.omit(smote_data_rf)
# Treinar modelo Random Forest
set.seed(123)
fit_rf4 <- train(Cancer ~ ., data = smote_data_rf, method = "rf",
trControl = ctrl1,
tuneGrid = expand.grid(mtry = 2), # Ajuste do número de variáveis
metric = "ROC",
ntree = 100) # Número de árvores ajustado
# Avaliação no teste
pred_rf <- predict(fit_rf4, testData)
conf_matrix <- confusionMatrix(pred_rf, testData$Cancer)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive:
table:
Astrocytoma Oligoastrocytoma Oligodendroglioma Astrocytoma 44 16 13 Oligoastrocytoma 46 29 25 Oligodendroglioma 1 17 48 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.5063 0.2694 0.4411 0.5713 0.3808 AccuracyPValue McnemarPValue 5.331e-05 8.151e-06 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Class: Astrocytoma 0.4835 0.8041 0.6027 Class: Oligoastrocytoma 0.4677 0.5989 0.29 Class: Oligodendroglioma 0.5581 0.8824 0.7273 Table continues below Neg Pred Value Precision Recall F1 Class: Astrocytoma 0.7169 0.6027 0.4835 0.5366 Class: Oligoastrocytoma 0.7626 0.29 0.4677 0.358 Class: Oligodendroglioma 0.7803 0.7273 0.5581 0.6316 Table continues below Prevalence Detection Rate Class: Astrocytoma 0.3808 0.1841 Class: Oligoastrocytoma 0.2594 0.1213 Class: Oligodendroglioma 0.3598 0.2008 Detection Prevalence Balanced Accuracy Class: Astrocytoma 0.3054 0.6438 Class: Oligoastrocytoma 0.4184 0.5333 Class: Oligodendroglioma 0.2762 0.7202 mode: sens_spec
dots:
Matriz de Confusão:
- Astrocytoma: 44 corretos, 16 classificados como Oligoastrocytoma, 13 como Oligodendroglioma.
- Oligoastrocytoma: 29 corretos, 46 classificados como Astrocytoma, 25 como Oligodendroglioma.
- Oligodendroglioma: 48 corretos, 17 classificados como Oligoastrocytoma, 1 como Astrocytoma.
Astrocytoma e Oligodendroglioma são razoavelmente bem classificados. Oligoastrocytoma tem o pior desempenho, frequentemente confundida com Astrocytoma. A distribuição de erros mostra alguma capacidade de distinção, mas com sobreposição nas características entre as classes.
Métricas de Desempenho:
- Accuracy: 0.5063 → O modelo classifica corretamente 50.63% dos exemplos, o que pode ser melhorado.
- Kappa: 0.2694 → Indica concordância moderada entre previsões e valores reais.
- Mcnemar’s Test P-Value: 8.151e-06 → P-valor baixo, indicando diferenças estatisticamente significativas entre algumas classes.
Especificidade: O modelo consegue distinguir melhor Oligodendroglioma das demais classes. - Astrocytoma: 0.8041 - Oligoastrocytoma: 0.5989 - Oligodendroglioma: 0.8824
Balanced Accuracy: Astrocytoma e Oligodendroglioma têm desempenho razoável, sendo que Oligoastrocytoma tem a pior acurácia balanceada, indicando que o modelo tem dificuldades em prever essa classe corretamente. - Astrocytoma: 0.6438 - Oligoastrocytoma: 0.5333 - Oligodendroglioma: 0.7202
# Tabela confusion matrix
conf_table <- conf_matrix$table
# Inicializa data frame para armazenar metricas
class_metrics <- data.frame(Class = levels(testData$Cancer),
Precision = numeric(length(levels(testData$Cancer))),
Recall = numeric(length(levels(testData$Cancer))),
F1_Score = numeric(length(levels(testData$Cancer))))
# Ciclo sobre cada classe para calcular precisão, sensibilidade e F1 score
for (class in levels(testData$Cancer)) {
# Verdadeiros Positivos (VP): Previsões corretas da classe
TP <- conf_table[class, class]
# Falsos Positivos (FP): Previsões da classe, mas pertencem a outra
FP <- sum(conf_table[, class]) - TP
# Falsos Negativos (FN): Pertencem à classe, mas foram previstos como outra
FN <- sum(conf_table[class, ]) - TP
# Precisão: VP / (VP + FP)
precision <- ifelse((TP + FP) == 0, NA, TP / (TP + FP))
# Sensibilidade (Recall): VP / (VP + FN)
recall <- ifelse((TP + FN) == 0, NA, TP / (TP + FN))
# F1 Score: 2 * (Precisão * Sensibilidade) / (Precisão + Sensibilidade)
f1 <- ifelse(is.na(precision) | is.na(recall) | (precision + recall) == 0, NA, 2 * (precision * recall) / (precision + recall))
# Guardar os resultados
class_metrics[class_metrics$Class == class, ] <- c(class, precision, recall, f1)
}
# Converter colunas numéricas
class_metrics[, 2:4] <- lapply(class_metrics[, 2:4], as.numeric)
# Round
class_metrics_rounded <- class_metrics
class_metrics_rounded[, 2:4] <- round(class_metrics_rounded[, 2:4], 3)
# Tabela
knitr::kable(class_metrics_rounded, caption = "Per-Class Performance Metrics (Precision, Recall, F1 Score)")| Class | Precision | Recall | F1_Score |
|---|---|---|---|
| Astrocytoma | 0.484 | 0.603 | 0.537 |
| Oligoastrocytoma | 0.468 | 0.290 | 0.358 |
| Oligodendroglioma | 0.558 | 0.727 | 0.632 |
Recall:
- Astrocytoma: 0.4835
- Oligoastrocytoma: 0.4677
- Oligodendroglioma: 0.5581
Precision:
- Astrocytoma: 0.6027
- Oligoastrocytoma: 0.29
- Oligodendroglioma: 0.7273
F1 Score:
- Astrocytoma: 0.5366
- Oligoastrocytoma: 0.358
- Oligodendroglioma: 0.6316
# Calcular AUC e Curva ROC
roc_curve <- roc(testData$Cancer, as.numeric(pred_rf))
auc_score <- auc(roc_curve)
print(paste("AUC Score:", auc_score))[1] “AUC Score: 0.679457639135059”
Forma da Curva:
- A curva segue um padrão escalonado, indicando que o modelo está a tomar decisões em diferentes limiares de classificação.
Comparação com um Classificador Aleatório:
A linha diagonal representa um classificador aleatório (AUC = 0.5).Como a curva do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score: 0.6794 → Tem uma boa capacidade de distinguir entre as classes, especialmente Astrocytoma e Oligodendroglioma. No entanto, o modelo não é considerado útil (só acima de 0.75).
# Ver importância dos genes
importance <- varImp(fit_rf4)
# Extração variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
126432 |
126432 |
100.00000 |
5906 |
5906 |
97.13694 |
253558 |
253558 |
88.28653 |
23421 |
23421 |
85.67398 |
3099 |
3099 |
82.24723 |
84902 |
84902 |
80.48595 |
339105 |
339105 |
78.50936 |
112616 |
112616 |
77.07250 |
387923 |
387923 |
75.00653 |
63940 |
63940 |
74.42168 |
4116 |
4116 |
72.35777 |
27439 |
27439 |
72.31948 |
51411 |
51411 |
71.14274 |
57657 |
57657 |
69.72485 |
55611 |
55611 |
66.72946 |
1489 |
1489 |
66.63683 |
2898 |
2898 |
66.11533 |
8453 |
8453 |
65.89152 |
114822 |
114822 |
65.68206 |
4848 |
4848 |
65.66334 |
Importância dos Genes:
- Genes mais relevantes: 126432 com um peso de 100 e Gene 5906 com peso 97.
- Os primeiros 5 genes têm valores relativamente altos (acima de 80), sugerindo que são os principais responsáveis pela decisão do modelo.
- O modelo pode estar dependente de poucos genes-chave, enquanto outros têm impacto reduzido. Isso pode indicar que algumas características podem ser removidas sem afetar o desempenho.
# Visualizar distribuição de probabilidades
prob_rf <- predict(fit_rf4, testData, type = "prob")
boxplot(prob_rf[,1] ~ testData$Cancer, col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Astrocytoma:
- Mediana alta, indicando previsões confiantes e sem outliers, pelo que o modelo é estável e com alta confiança para esta classe.
Oligoastrocytoma
- Mediana mais baixa do que Astrocytoma e sem outliers, pelo que há menor confiança média nas previsões para esta classe.
Oligodendroglioma
- Mediana mais baixa do que Oligoastrocytoma com baixa variablidade e outlier acima do 0.5. O modelo é de baixa confiância ao prever esta classe.
11.2 Support Vector Machine
## B. Support Vector Machine com SMOTE ---- SURVIVAL
# Preparar dados
common_samples <- intersect(colnames(expr_all), names(f_survival))
logCPM_sub <- expr_all[, common_samples]
labels <- f_survival[common_samples]
# Criar data frame com amostras como linhas e genes como colunas
data_model <- as.data.frame(t(logCPM_sub))
data_model$SurvivalStatus <- factor(labels)
levels(data_model$SurvivalStatus) <- make.names(levels(data_model$SurvivalStatus)) # tornar nomes válidos
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$SurvivalStatus, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Aplicar SMOTE
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
smote_result <- SMOTE(trainData_numeric, target = trainData$SurvivalStatus, K = 3, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "SurvivalStatus"
smote_data$SurvivalStatus <- factor(smote_data$SurvivalStatus,
levels = levels(trainData$SurvivalStatus))
# Treinar modelo SVM
set.seed(123)
fit_svm1 <- train(SurvivalStatus ~ ., data = smote_data, method = "svmRadial",
trControl = ctrl,
preProcess = c("center", "scale"),
tuneLength=5,
metric = "ROC")
# Previsões no conjunto de teste
# Garantir que os níveis são iguais
testData$SurvivalStatus <- factor(testData$SurvivalStatus,
levels = levels(fit_svm1$trainingData$.outcome))
pred_class <- predict(fit_svm1, testData)
pred_class <- factor(pred_class, levels = levels(testData$SurvivalStatus))
# Matriz de confusão
confusion <- confusionMatrix(pred_class, testData$SurvivalStatus)
pander::pander(confusion, caption = "Confusion Matrix Statistics (per class)")positive: X0.ALIVE.OR.DEAD.TUMOR.FREE
table:
Table continues below X0.ALIVE.OR.DEAD.TUMOR.FREE X0.ALIVE.OR.DEAD.TUMOR.FREE 177 X1.DEAD.WITH.TUMOR 8 X1.DEAD.WITH.TUMOR X0.ALIVE.OR.DEAD.TUMOR.FREE 43 X1.DEAD.WITH.TUMOR 12 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.7875 0.2253 0.7303 0.8375 0.7708 AccuracyPValue McnemarPValue 0.2989 1.927e-06 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0.9568 0.2182 0.8045 0.6 0.8045 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0.9568 0.8741 0.7708 0.7375 0.9167 Balanced Accuracy 0.5875 mode: sens_spec
dots:
Matriz de Confusão:
- 177 previsões corretas para “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 12 previsões corretas para “X1.DEAD.WITH.TUMOR”.
- 8 erros ao classificar “X1.DEAD.WITH.TUMOR” como “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 43 erros ao classificar “XO.ALIVE.OR.DEAD.TUMOR.FREE” como “X1.DEAD.WITH.TUMOR”
Modelo tem boa capacidade de prever a classe “XO.ALIVE.OR.DEAD.TUMOR.FREE”, mas mostra dificuldades consideráveis em identificar corretamente “X1.DEAD.WITH.TUMOR”.
Métricas de Desempenho:
- Accuracy: 0.7875 → Classifica corretamente na maioria dos casos.
- Kappa: 0.2253 → Indica um acordo fraco a moderado entre as previsões e os valores reais.
- Especificidade: 0.2182 → O modelo tem dificultade em acertar os casos negativos.
- Balanced Accuracy: 0.5875 → O modelo é melhor do que o aleatório, mas ainda pode ser melhorado.
# Precision, Recall, F1
positive_class <- "X0.ALIVE.OR.DEAD.TUMOR.FREE"
precision <- posPredValue(pred_class, testData$SurvivalStatus, positive = positive_class)
recall <- sensitivity(pred_class, testData$SurvivalStatus, positive = positive_class)
f1 <- 2 * (precision * recall) / (precision + recall)
cat("Precision:", precision, "\n")Precision: 0.8045455
Recall: 0.9567568
F1 Score: 0.874
- Recall: 0.9568 → Identifica 95.68% dos casos positivos corretamente.
- Precisão: 0.8045 → Quando prevê “XO.ALIVE.OR.DEAD.TUMOR.FREE”, está correto 80.45% das vezes.
- F1 Score: 0.8741 → Indica um bom equilíbrio entre precisão e recall, mas reflete sobretudo a classe dominante.
# Prever probabilidades
pred_prob <- predict(fit_svm1, testData, type = "prob")
# AUC
correct_col <- grep("ALIVE.OR.DEAD.TUMOR.FREE", colnames(pred_prob), value = TRUE)
roc_curve <- roc(response = testData$SurvivalStatus, predictor = pred_prob[[correct_col]])
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.7808354
Forma da Curva:
- A curva ROC está num padrão escalonado, indicando que o modelo está a tomar decisões em diferentes limiares de classificação.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5).
- Como a curva ROC do modelo está bem acima dessa linha, isso confirma que o modelo tem valor preditivo real.
AUC Score: - AUC: 0.7808 → O modelo tem capacidade razoável de distinguir entre as classes, mas não é perfeito (abaixo de 0.8).
# Variavel importance
importance <- varImp(fit_svm1, scale = FALSE)
# importance em data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df$Overall <- rowMeans(importance_df[, c("X0.ALIVE.OR.DEAD.TUMOR.FREE", "X1.DEAD.WITH.TUMOR")])
# Sort por importance
importance_df <- importance_df[order(-importance_df$Overall), ]
# Seleciona top 20
top_20 <- head(importance_df, 20)
# tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
| 358 | 358 | 0.7792766 |
| 63877 | 63877 | 0.7680026 |
| 5738 | 5738 | 0.7642880 |
| 284805 | 284805 | 0.7627566 |
| 3908 | 3908 | 0.7621701 |
| 56852 | 56852 | 0.7615836 |
| 128977 | 128977 | 0.7604431 |
| 8786 | 8786 | 0.7594982 |
| 26524 | 26524 | 0.7583904 |
| 9180 | 9180 | 0.7581297 |
| 7764 | 7764 | 0.7577061 |
| 55970 | 55970 | 0.7571848 |
| 64778 | 64778 | 0.7546758 |
| 135228 | 135228 | 0.7521342 |
| 7060 | 7060 | 0.7521017 |
| 57715 | 57715 | 0.7520691 |
| 3485 | 3485 | 0.7520365 |
| 64114 | 64114 | 0.7520365 |
| 1116 | 1116 | 0.7514174 |
| 202333 | 202333 | 0.7497882 |
Importância dos Genes:
- Gene mais relevante: 358 com um peso de 0.07793. Este gene tem um impacto significativo na classificação do modelo.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.07), sugerindo que são os principais responsáveis pela decisão do modelo
# Boxplot da distribuição de probabilidades
df_plot <- data.frame(Prob = pred_prob[[correct_col]], Status = testData$SurvivalStatus)
boxplot(Prob ~ Status, data = df_plot,
col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe “X0.ALIVE.OR.DEAD.TUMOR.FREE” (vermelho):
- A maioria das previsões tem probabilidades próximas de 1, indicando alta confiança do modelo nesta classe.
- Existem alguns outliers com probabilidades mais baixas (~0.9), sugerindo casos com menor certeza, possivelmente mais difíceis de classificar corretamente.
Classe “X1.DEAD.WITH.TUMOR” (azul):
- A distribuição é mais dispersa, com valores variando de 0.5 a 0.9.
- A mediana está próxima de 0.8, indicando que o modelo não tem tanta certeza ao prever esta classe. Pode sugerir dificuldade na separação entre as classes, levando a previsões menos confiáveis.
## B. Support Vector Machine com SMOTE ---- SEX
# Preparar dados
common_samples <- intersect(colnames(expr_all), names(f_sex))
logCPM_sub <- expr_all[, common_samples]
labels <- f_sex[common_samples]
# Criar data frame com amostras como linhas e genes como colunas
data_model <- as.data.frame(t(logCPM_sub))
data_model$Sex <- factor(labels)
# Garantir que os levels são válidos para uso com class probabilities
levels(data_model$Sex) <- make.names(levels(data_model$Sex))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Sex, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Corrigir aplicação do SMOTE
# Remover colunas não numéricas
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
# Encontrar o tamanho mínimo da classe
min_class_size <- min(table(trainData$Sex))
# Ajustar K para evitar erro de SMOTE
K_value <- max(1, min_class_size - 1)
# Aplicar SMOTE corretamente (usando factor numérico)
smote_result <- SMOTE(trainData_numeric, target = trainData$Sex, K = K_value, dup_size = 2)
smote_data <- smote_result$data
# Ajustar nomes das colunas
colnames(smote_data)[ncol(smote_data)] <- "Sex"
smote_data$Sex <- factor(smote_data$Sex, levels = levels(trainData$Sex))
# Remover valores ausentes
smote_data <- na.omit(smote_data)
# Treinar modelo SVM com validação cruzada e métrica ROC
set.seed(123)
fit_svm2 <- train(Sex ~ ., data = smote_data, method = "svmRadial",
trControl = ctrl,
preProcess = c("center", "scale"),
tuneLength=5,
metric = "ROC")
# Definir a classe positiva para avaliação
positive_class <- "Female"
# Previsão do modelo no conjunto de teste
pred_class <- predict(fit_svm2, testData)
# Criar matriz de confusão para avaliar o desempenho do modelo
confusion <- confusionMatrix(pred_class, testData$Sex)
pander::pander(confusion, caption = "Estatísticas da Matriz de Confusão (por classe)")positive: Female
table:
Female Male Female 36 16 Male 70 118 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.6417 0.2326 0.5775 0.7023 0.5583 AccuracyPValue McnemarPValue 0.005329 1.096e-08 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0.3396 0.8806 0.6923 0.6277 0.6923 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0.3396 0.4557 0.4417 0.15 0.2167 Balanced Accuracy 0.6101 mode: sens_spec
dots:
# Garantir que testData$Sex é um fator
testData$Sex <- factor(testData$Sex)
# Garantir que as previsões têm os mesmos níveis que os dados reais
pred_svm <- factor(pred_class, levels = levels(testData$Sex)) Matriz de confusão:
- 36 previsões corretas para “Female” e 118 previsões corretas para “Male”.
- 16 erros ao classificar “Female” como “Male”.
- 70 erros ao classificar “Male” como “Female”.
Modelo apresenta bom desempenho na classe “Male”, mas desempenho fraco na classe “Female”. A aplicação de SMOTE não foi suficiente para corrigir o desequílibrio nem para melhorar de forma significativa a sensibilidade para a classe minoritária.
Métricas de Desempenho:
- Accuracy: 0.6417 → o modelo acerta aproximadamente dois terços dos casos.
- Kappa: 0.2326 → desempenho ligeiramente melhor do modelo aleatório.
- Especificidade: 0.8806 → o modelo reconhece bem a classe “Male”, com poucos falsos positivos.
# Calcular Precisão, Recall e F1 Score
precision <- posPredValue(pred_svm, testData$Sex, positive = positive_class)
recall <- sensitivity(pred_svm, testData$Sex, positive = positive_class)
# Exibir os valores calculados
cat("Precisão:", precision, "\n")Precisão: 0.6923077
Recall: 0.3396226
# Calcular o F1 Score apenas se houver valores válidos de precisão e recall
if (!is.na(precision) && !is.na(recall) && (precision + recall) > 0) {
f1 <- 2 * (precision * recall) / (precision + recall)
cat("F1 Score:", round(f1, 3), "\n")
} else {
cat("O F1 Score não pode ser calculado devido à ausência de precisão ou recall.\n")
}F1 Score: 0.456
# Previsão das probabilidades para cada classe
pred_prob <- predict(fit_svm2, testData, type = "prob")
# Garantir que testData$Sex tem os mesmos níveis que smote_data$Sex
testData$Sex <- factor(testData$Sex, levels = levels(smote_data$Sex))
# Encontrar dinamicamente a coluna correta de probabilidades
sex_levels <- levels(testData$Sex)
correct_col <- grep(paste(sex_levels, collapse = "|"), colnames(pred_prob), value = TRUE)
# Verificar se a coluna correta foi encontrada
if (length(correct_col) == 0) {
stop(paste("Erro: Coluna de probabilidade para a classificação de Sex não encontrada. Nomes disponíveis:",
paste(colnames(pred_prob), collapse = ", ")))
}
# Garantir que apenas uma coluna é selecionada
if (length(correct_col) > 1) {
correct_col <- correct_col[1] # Selecionar a primeira correspondência
}
# Garantir que o preditor e a resposta têm o mesmo comprimento
if (length(testData$Sex) != nrow(pred_prob)) {
stop("Erro: Incompatibilidade de tamanhos entre testData$Sex e as probabilidades previstas.")
}
# Converter preditor para numérico
numeric_pred <- as.numeric(pred_prob[[correct_col]]) - Precisão: 0.6923 → Quando o modelo prevê “Female”, está correto 69.23% das vezes.
- Recall: 0.3396 → Identifica 33.96% dos casos “Female” corretamente.
- F1 Score: 0.4557 → Indica baixo equilíbrio entre precisão e recall para esta classe.
# Calcular a curva ROC
roc_curve <- roc(response = testData$Sex, predictor = numeric_pred)
# Calcular a AUC (Área sob a curva ROC)
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.7230358
Forma da Curva:
- A curva está num padrão escalonado, indicando que o modelo está a tomar decisões com base em múltiplos limiares de classificação.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva do modelo está acima da diagonal, o modelo tem desempenho superior ao acaso.
AUC Score:
- AUC: 0.723 → O modelo tem capacidade moderada de distinguir entre as classes “Female” e “Male”.
# Calcular a importância das variáveis (genes)
importance <- varImp(fit_svm2, scale = FALSE)
# Converter importância para um data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
# Criar uma métrica geral de importância combinando as classes "Female" e "Male"
importance_df$Overall <- rowMeans(importance_df[, c("Female", "Male")])
# Ordenar por importância
importance_df <- importance_df[order(-importance_df$Overall), ]
# Selecionar as 20 características mais importantes
top_20 <- head(importance_df, 20)
# Exibir tabela com os genes mais importantes
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Genes Mais Importantes")| Gene | Overall | |
|---|---|---|
| 22829 | 22829 | 0.9980239 |
| 554203 | 554203 | 0.8700423 |
| 8242 | 8242 | 0.8629051 |
| 7403 | 7403 | 0.8539313 |
| 7543 | 7543 | 0.8349607 |
| 8226 | 8226 | 0.8246152 |
| 1654 | 1654 | 0.7763519 |
| 7317 | 7317 | 0.7744455 |
| 94056 | 94056 | 0.7697726 |
| 1968 | 1968 | 0.7583577 |
| 1964 | 1964 | 0.7558469 |
| 9189 | 9189 | 0.7374576 |
| 8623 | 8623 | 0.7246710 |
| 412 | 412 | 0.7079323 |
| 55787 | 55787 | 0.6823127 |
| 8239 | 8239 | 0.6817315 |
| 80230 | 80230 | 0.6774771 |
| 207063 | 207063 | 0.6727112 |
| 3563 | 3563 | 0.6648533 |
| 4267 | 4267 | 0.6634817 |
# Plot da importância das características
plot(importance, top = 20, main = "Top 20 Características Mais Importantes")
Importância dos Genes:
- Gene mais relevante: 22829 com um peso de 0.0998, o qual domina completamente a decisão do modelo.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.08), sugerindo que são os principais responsáveis pela decisão do modelo.
# Criar um data frame para o boxplot
df_plot <- data.frame(Prob = pred_prob[, correct_col], Status = testData$Sex)
# Gerar o boxplot para visualizar a distribuição das probabilidades previstas
boxplot(pred_prob[, correct_col] ~ testData$Sex,
col = palette_3,
main = "Distribuição das Probabilidades por Classe", ylab = "Probabilidade Prevista")
Boxplot - Distribuição das Probabilidades:
Classe “Female” (vermelho):
- Possui uma mediana perto dos 0.3, indicando previsões moderadamente confiantes.
- Dispersão alta com o Q3 perto dos 0.6 e o Q1 perto dos 0.1, isto sugere que o modelo tem dificuldade em ser consistente ao prever esta classe.
Classe “Male” (azul):
- Possui uma mediana baixa perto dos 0.05, indicando previsões mais consistentes para prever esta classe.
- Poucos outliers acima dos 0.6, sugerindo robustez do modelo.
A sobreposição nas distribuições de probabilidade entre as duas classes mostra que o modelo tem dificuldade em separá-las com confiança. Isto confirma o desempenho desigual observado nas métricas e matriz de confusão, particularmente a fraca performance na classe “Female”.
# --- B. Support Vector Machine com SMOTE ---- ANCESTRALIDADE
# Preparar dados
clean_data_for_2part <- clini_data_no_na
# Definir variável resposta
f_ancestralidade <- as.factor(clini_data_no_na$ancestry)
names(f_ancestralidade) <- clini_data_no_na$sample_ID
# Alinhar dados de expressão com amostras com dados clínicos
common_samples <- intersect(colnames(expr_all), names(f_ancestralidade))
logCPM_sub <- expr_all[, common_samples]
labels <- f_ancestralidade[common_samples]
# Criar data frame com amostras como linhas e genes como colunas
data_model <- as.data.frame(t(logCPM_sub))
data_model$Ancestralidade <- factor(labels)
# Garantir que os levels são válidos para uso com class probabilities
levels(data_model$Ancestralidade) <- make.names(levels(data_model$Ancestralidade))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Ancestralidade, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# --- FILTRAR CLASSES RARAS (< 5 amostras) ---
min_samples <- 5
table_counts <- table(trainData$Ancestralidade)
keep_levels <- names(table_counts[table_counts >= min_samples])
trainData_filtered <- trainData[trainData$Ancestralidade %in% keep_levels, ]
trainData_filtered$Ancestralidade <- droplevels(trainData_filtered$Ancestralidade)
# Atualizar levels no conjunto de teste para corresponder
testData <- testData[testData$Ancestralidade %in% keep_levels, ]
testData$Ancestralidade <- droplevels(testData$Ancestralidade)
# Preparar dados numéricos para SMOTE
trainData_numeric <- trainData_filtered[, sapply(trainData_filtered, is.numeric)]
# Encontrar novo tamanho mínimo da classe
min_class_size <- min(table(trainData_filtered$Ancestralidade))
K_value <- max(1, min_class_size - 1)
# Aplicar SMOTE
smote_result <- SMOTE(trainData_numeric, target = trainData_filtered$Ancestralidade, K = K_value, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "Ancestralidade"
smote_data$Ancestralidade <- factor(smote_data$Ancestralidade, levels = levels(trainData_filtered$Ancestralidade))
# Remover NAs
smote_data <- na.omit(smote_data)
# Treinar modelo SVM
set.seed(123)
fit_svm3 <- train(Ancestralidade ~ ., data = smote_data, method = "svmRadial",
trControl = ctrl1,
preProcess = c("center", "scale"),
tuneLength=5,
metric = "Accuracy")
# Previsão do modelo no conjunto de teste
pred_class <- predict(fit_svm3, testData)
# Criar matriz de confusão
conf_matrix <- confusionMatrix(pred_class, testData$Ancestralidade)
pander::pander(conf_matrix, caption = "Estatísticas da Matriz de Confusão (por classe)")positive: AFR
table:
AFR EUR AFR 0 0 EUR 7 218 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.9689 0 0.937 0.9874 0.9689 AccuracyPValue McnemarPValue 0.5987 0.02334 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0 1 NA 0.9689 NA Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0 NA 0.03111 0 0 Balanced Accuracy 0.5 mode: sens_spec
dots:
Matriz de Confusão:
- 218 previsões corretas para “EUR”, mas teve dificuldades nas outras classes.
- Confirma que o modelo está enviesado para prever EUR corretamente, mas falha nas outras classes
Métricas de Desempenho:
- Accuracy: 0.9689 → Classifica corretamente na maioria dos casos para a classe dominante.
- AccuracyLower: 0.937 e AccuracyUpper: 0.9874 → o intervalo de confiança de accuracy mostra estabilidade na mesma tendência enviesada
- AccuracyNull: 0.9689 → o modelo pode estar apenas a prever a classe dominante.
- Kappa: 0 → indica que o modelo não está a prever melhor do que um classificador aleatório.
- Especificidade: 1 → como o modelo não identifica corretamente as classes negativas, parece haver confusão com a sensibilidade ou precisão.
- Balanced Accuracy: 0.5 → sugere um desempenho aleatório quando se consideram todas as classes iguais
- p-value da accuracy (0.5987) indica que não há diferença estatística significativa entre o modelo e um classificador aleatório.
- McNemarPValue: 0.02334 → possível desequilíbrio nos erros, isto é, o modelo trata uma classe de forma muito diferente das outras.
# Calcular Precisão, Recall e F1 Score por classe
precision_vec <- diag(conf_matrix$table) / colSums(conf_matrix$table)
recall_vec <- diag(conf_matrix$table) / rowSums(conf_matrix$table)
# Calcular F1-score para cada classe
f1_vec <- ifelse((precision_vec + recall_vec) > 0,
2 * (precision_vec * recall_vec) / (precision_vec + recall_vec),
NA)
# Criar um data frame com os valores de F1-score
f1_scores <- data.frame(
Classe = rownames(conf_matrix$table),
Precisão = round(precision_vec, 3),
Recall = round(recall_vec, 3),
F1 = round(f1_vec, 3)
)
cat("F1, Recall e Precisão por classe:\n")F1, Recall e Precisão por classe:
| Classe | Precisão | Recall | F1 | |
|---|---|---|---|---|
| AFR | AFR | 0 | NaN | NA |
| EUR | EUR | 1 | 0.969 | 0.984 |
# Previsão de probabilidades (para AUC e boxplot)
pred_prob <- predict(fit_svm3, testData, type = "prob")
# Garantir que testData$Ancestralidade tem os mesmos níveis que pred_prob
testData$Ancestralidade <- factor(testData$Ancestralidade, levels = colnames(pred_prob))- Precisão: NA → Prevê apenas a classe “EUR” com alta confiança.
- Recall: 0 → Identifica bem apenas a classe “EUR”.
- F1 Score: NA → Prevê apenas a classe “EUR” com alta confiança.
# Criar listas para armazenar curvas ROC e valores de AUC
roc_list <- list()
auc_scores <- c()
classes <- colnames(pred_prob)
# Calcular curvas ROC para cada classe (um-contra-todos)
for (class_name in classes) {
# Criar rótulos binários: 1 para a classe atual, 0 para as restantes
binary_response <- factor(ifelse(testData$Ancestralidade == class_name, class_name, paste0("not_", class_name)))
# Calcular curva ROC
roc_curve <- roc(response = binary_response,
predictor = pred_prob[, class_name],
levels = c(paste0("not_", class_name), class_name))
# Guardar valores de AUC e curva ROC
auc_score <- auc(roc_curve)
auc_scores <- c(auc_scores, auc_score)
roc_list[[class_name]] <- roc_curve
cat(paste("AUC para", class_name, "vs todos:", round(auc_score, 3), "\n"))
}AUC para AFR vs todos: 0.649 AUC para EUR vs todos: 0.649
# Plot a primeira curva ROC
plot(roc_list[[1]], col = palette_3[1], main = "Curvas ROC Multi-classe para SVM (com SMOTE)")
Forma da Curva:
- A curva ROC tem uma forma degraudo e irregular, permanecendo perto da linha diagonal, o que indica baixa capacidade de separação entre classes.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva está um pouco acima da diagonal, o desempenho é ligeiramente melhor que o acaso, mas sem ganhos significativos.
AUC Score:
- AUC: 0.649 → o modelo tem um desempenho fraco na sua capacidade discriminativa, mas é melhor que aleatório.
# Calcular a importância das variáveis (genes)
importance <- varImp(fit_svm3, scale = FALSE)
# Converter importância para um data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
# Criar uma métrica geral de importância
importance_df$Overall <- rowMeans(importance_df[, c("AFR", "EUR")])
# Ordenar por importância
importance_df <- importance_df[order(-importance_df$Overall), ]
# Selecionar as 20 características mais importantes
top_20 <- head(importance_df, 20)
# Exibir tabela com os genes mais importantes
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Genes Mais Importantes")| Gene | Overall | |
|---|---|---|
| 317749 | 317749 | 0.9575688 |
| 10901 | 10901 | 0.9248853 |
| 84221 | 84221 | 0.9099771 |
| 100132288 | 100132288 | 0.8965979 |
| 54986 | 54986 | 0.8948777 |
| 79446 | 79446 | 0.8814985 |
| 9652 | 9652 | 0.8700306 |
| 93100 | 93100 | 0.8635321 |
| 6884 | 6884 | 0.8627676 |
| 3920 | 3920 | 0.8604740 |
| 2979 | 2979 | 0.8576070 |
| 83787 | 83787 | 0.8564602 |
| 115752 | 115752 | 0.8562691 |
| 23061 | 23061 | 0.8518731 |
| 2318 | 2318 | 0.8511086 |
| 606495 | 606495 | 0.8503440 |
| 64784 | 64784 | 0.8488150 |
| 64397 | 64397 | 0.8478593 |
| 2298 | 2298 | 0.8476682 |
| 162394 | 162394 | 0.8470948 |
# Plotar a importância das características
plot(importance, top = 20, main = "Top 20 Características Mais Importantes")
Importância dos Genes:
- Genes mais relevantes: 317749 com um peso de 0.9576 e Gene 10901 com peso 0.9249.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.09), sugerindo que são os principais responsáveis pela decisão do modelo.
### Correção do Boxplot da Distribuição de Probabilidades
# Encontrar dinamicamente a coluna correta de probabilidades
correct_col1 <- grep(paste(levels(testData$Ancestralidade), collapse = "|"), colnames(pred_prob), value = TRUE)
# Garantir que apenas uma coluna é selecionada
if (length(correct_col1) > 1) {
correct_col1 <- correct_col1[1] # Selecionar a primeira correspondência
}
# Verificar se a coluna correta foi encontrada
if (!(correct_col1 %in% colnames(pred_prob))) {
stop("Erro: Coluna de probabilidade para a classificação de Ancestralidade não encontrada.")
}
# Converter preditor para numérico
numeric_pred <- as.numeric(pred_prob[[correct_col1]])
# Criar data frame para o boxplot
df_plot <- data.frame(Prob = numeric_pred, Status = testData$Ancestralidade)
# Gerar o boxplot
boxplot(df_plot$Prob ~ df_plot$Status,
col = palette_3,
main = "Distribuição das Probabilidades por Classe", ylab = "Probabilidade Prevista")
Boxplot - Distribuição das Probabilidades:
Classe AFR:
- A mediana da probabilidade predita está próxima de 0.01.
- Existe uma baixa dispersão e não há outliers, sugerindo que o modelo prevê esta classe de forma repetitiva, mas com baixa confiança, indicando a sua incapacidade de discriminação real.
Classe EUR:
- A mediana também é próxima de 0.01, e possui outliers acima de 0.02, indicando alguns casos onde o modelo atribui maior probabilidade.
- A presença de outliers em EUR mostra tendência a superestimar a classe dominante mesmo que discretamente.
## B. Support Vector Machine com SMOTE ---- CANCER
# Preparar dados
common_samples <- intersect(colnames(expr_all), names(f_cancer))
logCPM_sub <- expr_all[, common_samples]
labels <- f_cancer[common_samples]
# Criar data frame com amostras como linhas e genes como colunas
data_model <- as.data.frame(t(logCPM_sub))
data_model$Cancer <- factor(labels)
# Garantir que os levels são válidos para uso com class probabilities
levels(data_model$Cancer) <- make.names(levels(data_model$Cancer))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Cancer, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Remover colunas não numéricas
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
# Encontrar o tamanho mínimo da classe
min_class_size <- min(table(trainData$Cancer))
# Ajustar K para evitar erro de SMOTE
K_value <- max(1, min_class_size - 1) # Garante que K seja pelo menos 1
# Aplicar SMOTE corretamente (usando factor numérico)
smote_result <- SMOTE(trainData_numeric, target = trainData$Cancer, K = K_value, dup_size = 2)
smote_data <- smote_result$data
# Ajustar nomes das colunas
colnames(smote_data)[ncol(smote_data)] <- "Cancer"
smote_data$Cancer <- factor(smote_data$Cancer, levels = levels(trainData$Cancer))
# Remover valores ausentes
smote_data <- na.omit(smote_data)
# Treinar modelo SVM com validação cruzada e métrica ROC
set.seed(123)
fit_svm4 <- train(Cancer ~ ., data = smote_data, method = "svmRadial",
trControl = ctrl1, # Use multiClassSummary
preProcess = c("center", "scale"),
tuneLength=5,
metric = "Accuracy") # Use Accuracy for multi-class problems
# Avaliação no conjunto de teste
pred_class <- predict(fit_svm4, testData)
confusion <- confusionMatrix(pred_class, testData$Cancer)
pander::pander(confusion, caption = "Confusion Matrix Statistics (per class)")positive:
table:
Astrocytoma Oligoastrocytoma Oligodendroglioma Astrocytoma 69 31 21 Oligoastrocytoma 8 9 3 Oligodendroglioma 14 22 62 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.5858 0.3507 0.5205 0.6489 0.3808 AccuracyPValue McnemarPValue 1.138e-10 1.842e-06 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Class: Astrocytoma 0.7582 0.6486 0.5702 Class: Oligoastrocytoma 0.1452 0.9379 0.45 Class: Oligodendroglioma 0.7209 0.7647 0.6327 Table continues below Neg Pred Value Precision Recall F1 Class: Astrocytoma 0.8136 0.5702 0.7582 0.6509 Class: Oligoastrocytoma 0.758 0.45 0.1452 0.2195 Class: Oligodendroglioma 0.8298 0.6327 0.7209 0.6739 Table continues below Prevalence Detection Rate Class: Astrocytoma 0.3808 0.2887 Class: Oligoastrocytoma 0.2594 0.03766 Class: Oligodendroglioma 0.3598 0.2594 Detection Prevalence Balanced Accuracy Class: Astrocytoma 0.5063 0.7034 Class: Oligoastrocytoma 0.08368 0.5415 Class: Oligodendroglioma 0.41 0.7428 mode: sens_spec
dots:
Matriz de Confusão:
- Astrocytoma: 69 corretos, 31 classificados como Oligoastrocytoma, 21 como Oligodendroglioma.
- Oligoastrocytoma: 9 corretos, 8 classificados como Astrocytoma, 3 como Oligodendroglioma.
- Oligodendroglioma: 62 corretos, 22 classificados como Oligoastrocytoma, 14 como Astrocytoma. Astrocytoma e Oligodendroglioma são razoavelmente bem classificados. Oligoastrocytoma tem baixo desempenho, frequentemente confundida com outras. Isto sugere sobreposição nas características dos tumores e limites poucos claros entre as classes.
Métricas de Desempenho:
- Accuracy: 0.5858 → O modelo classifica corretamente 58.58% dos exemplos, o que pode ser melhorado.
- Kappa: 0.3507 → Indica concordância moderada entre previsões e valores reais.
- Mcnemar’s Test P-Value: 1.842e-06 → P-valor baixo, indicando diferenças estatisticamente significativas entre algumas classes.
Especificidade: O modelo consegue distinguir melhor Astrocytoma e Oligodendroglioma. - Astrocytoma: 0.6486 - Oligoastrocytoma: 0.9379 - Oligodendroglioma: 0.7647
Balanced Accuracy: Astrocytoma e Oligodendroglioma têm desempenho razoável, sendo que Oligoastrocytoma tem a pior acurácia balanceada, indicando que o modelo tem dificuldades em prever essa classe corretamente. - Astrocytoma: 0.7034 - Oligoastrocytoma: 0.5415 - Oligodendroglioma: 0.7428
# Extrair tabela confusion matrix
conf_table <- confusion$table
# Initializa data frame para armazenar métricas
class_metrics <- data.frame(Class = levels(testData$Cancer),
Precision = numeric(length(levels(testData$Cancer))),
Recall = numeric(length(levels(testData$Cancer))),
F1_Score = numeric(length(levels(testData$Cancer))))
# Ciclo sobre cada classe para calcular precisão, sensibilidade e F1 score
for (class in levels(testData$Cancer)) {
# Verdadeiros Positivos (VP): Previsões corretas da classe
TP <- conf_table[class, class]
# Falsos Positivos (FP): Previsões da classe, mas pertencem a outra
FP <- sum(conf_table[, class]) - TP
# Falsos Negativos (FN): Pertencem à classe, mas foram previstos como outra
FN <- sum(conf_table[class, ]) - TP
# Precisão: VP / (VP + FP)
precision <- ifelse((TP + FP) == 0, NA, TP / (TP + FP))
# Sensibilidade (Recall): VP / (VP + FN)
recall <- ifelse((TP + FN) == 0, NA, TP / (TP + FN))
# F1 Score: 2 * (Precisão * Sensibilidade) / (Precisão + Sensibilidade)
f1 <- ifelse(is.na(precision) | is.na(recall) | (precision + recall) == 0, NA, 2 * (precision * recall) / (precision + recall))
# Guardar os resultados
class_metrics[class_metrics$Class == class, ] <- c(class, precision, recall, f1)
}
# Converter colunas numéricas
class_metrics[, 2:4] <- lapply(class_metrics[, 2:4], as.numeric)
# Round
class_metrics_rounded <- class_metrics
class_metrics_rounded[, 2:4] <- round(class_metrics_rounded[, 2:4], 3)
# Tabela
knitr::kable(class_metrics_rounded, caption = "Per-Class Performance Metrics (Precision, Recall, F1 Score)")| Class | Precision | Recall | F1_Score |
|---|---|---|---|
| Astrocytoma | 0.758 | 0.570 | 0.651 |
| Oligoastrocytoma | 0.145 | 0.450 | 0.220 |
| Oligodendroglioma | 0.721 | 0.633 | 0.674 |
# Prever probabilidades
pred_prob <- predict(fit_svm4, testData, type = "prob")
# Garantir que testData$Cancer tem os mesmos níveis que smote_data$Cancer
testData$Cancer <- factor(testData$Cancer, levels = levels(smote_data$Cancer))
# Encontrar a coluna correta para probabilidades
cancer_levels <- levels(testData$Cancer)
correct_col <- grep(paste(cancer_levels, collapse = "|"), colnames(pred_prob), value = TRUE)
# Verificar se a coluna correta foi encontrada
if (length(correct_col) == 0) {
stop(paste("Erro: Coluna de probabilidade para classificação de Cancer não encontrada. Nomes disponíveis:",
paste(colnames(pred_prob), collapse = ", ")))
}
# Garantir que apenas uma coluna seja selecionada
if (length(correct_col) > 1) {
correct_col <- correct_col[1]
}
# Verificar se os tamanhos das variáveis são compatíveis
if (length(testData$Cancer) != nrow(pred_prob)) {
stop("Erro: Diferente número de elementos entre testData$Cancer e pred_prob.")
}Recall:
- Astrocytoma: 0.7582
- Oligoastrocytoma: 0.1452
- Oligodendroglioma: 0.7209
Precision:
- Astrocytoma: 0.5702
- Oligoastrocytoma: 0.45
- Oligodendroglioma: 0.6327
F1 Score: - Astrocytoma: 0.6509 - Oligoastrocytoma: 0.2195 - Oligodendroglioma: 0.6739
# Converter preditor para numérico
numeric_pred <- as.numeric(pred_prob[[correct_col]]) # Usa [[ ]] para extrair como vetor
# Calcular AUC e Curva ROC
roc_curve <- roc(response = testData$Cancer, predictor = numeric_pred)
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.7227933
Forma da Curva:
- A curva segue um padrão escalonado, indicando que o modelo está a tomar decisões em diferentes limiares de classificação.
Comparação com um Classificador Aleatório:
A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score: 0.72279 → Tem uma capacidade moderada de distinguir entre as classes, especialmente Astrocytoma e Oligodendroglioma.
# Calcular a importância das variáveis no modelo
importance <- varImp(fit_svm4, scale = FALSE)
# Converter a importância para um data frame
importance_df <- as.data.frame(importance$importance)
# Adicionar os nomes dos genes como uma coluna separada
importance_df$Gene <- rownames(importance_df)
# Criar uma métrica geral de importância combinando as classes "Atrocytoma", "Oligoastrocytoma" e "Oligodendroglioma"
importance_df$Overall <- rowMeans(importance_df[, c("Astrocytoma", "Oligoastrocytoma", "Oligodendroglioma")])
# Ordenar os genes por importância (do maior para o menor)
importance_df <- importance_df[order(-importance_df$Overall), ]
# Selecionar os 20 genes mais importantes
top_20 <- head(importance_df, 20)
# Exibir tabela com os genes mais importantes
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Genes Mais Importantes")| Gene | Overall | |
|---|---|---|
| 83931 | 83931 | 0.8397411 |
| 55616 | 55616 | 0.8385525 |
| 79971 | 79971 | 0.8375420 |
| 204 | 204 | 0.8346363 |
| 1266 | 1266 | 0.8292582 |
| 51075 | 51075 | 0.8290788 |
| 832 | 832 | 0.8264122 |
| 7431 | 7431 | 0.8254620 |
| 1114 | 1114 | 0.8241157 |
| 85440 | 85440 | 0.8227958 |
| 10487 | 10487 | 0.8223648 |
| 3065 | 3065 | 0.8220713 |
| 27095 | 27095 | 0.8218554 |
| 9957 | 9957 | 0.8203804 |
| 65979 | 65979 | 0.8202015 |
| 8915 | 8915 | 0.8195505 |
| 2787 | 2787 | 0.8194699 |
| 200081 | 200081 | 0.8183607 |
| 3767 | 3767 | 0.8182258 |
| 79084 | 79084 | 0.8177135 |
# Plotar a importância das características
plot(importance, top = 20, main = "Top 20 Características Mais Importantes")
Importância dos Genes:
- Genes mais relevantes: 83931 com um peso de 0.8397 e Gene 55616 com peso 0.8385.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.8), sugerindo que são os principais responsáveis pela decisão do modelo.
# Boxplot da distribuição de probabilidades
df_plot <- data.frame(Prob = pred_prob[[correct_col]], Status = testData$Cancer)
boxplot(pred_prob[[correct_col]] ~ testData$Cancer,
col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Astrocytoma:
- Mediana alta, indicando previsões confiantes e sem outliers, pelo que o modelo é estável e com alta confiança para esta classe.
Oligoastrocytoma:
- Mediana mais baixa do que Astrocytoma e sem outliers, pelo que há menor confiança média nas previsões para esta classe.
Oligodendroglioma:
- Mediana mais baixa do que Oligoastrocytoma com baixa variablidade e sem outliers. O modelo é de baixa confiância ao prever esta classe.
11.3 XGBoost
## C. XGBoost com SMOTE ---- SURVIVAL
# Preparar dados
clean_data_for_2part <- clini_data_no_na
# Definir variável resposta
f_survival <- as.factor(clini_data_no_na$surv_status)
names(f_survival) <- clini_data_no_na$sample_ID
# Alinhar dados de expressão com amostras com dados clínicos
common_samples <- intersect(colnames(expr_all), names(f_survival))
logCPM_sub <- expr_all[, common_samples]
labels <- f_survival[common_samples]
# Construir dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub)) # genes como colunas
data_model$SurvivalStatus <- factor(labels)
# Ajustar levels para uso com SMOTE e xgboost
levels(data_model$SurvivalStatus) <- make.names(levels(data_model$SurvivalStatus))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$SurvivalStatus, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# SMOTE
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
smote_result <- SMOTE(trainData_numeric, trainData$SurvivalStatus, K = 3, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "SurvivalStatus"
smote_data$SurvivalStatus <- factor(smote_data$SurvivalStatus,
levels = levels(trainData$SurvivalStatus))
# Treinar modelo XGBoost
set.seed(123)
capture.output({
fit_xgb1 <- train(SurvivalStatus ~ ., data = smote_data, method = "xgbTree",
trControl = ctrl,
metric ="ROC",
tuneLength=3)
}, file = "NUL")
# Avaliação no conjunto de teste
pred_class <- predict(fit_xgb1, testData)
conf_matrix <- confusionMatrix(pred_class, testData$SurvivalStatus)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive: X0.ALIVE.OR.DEAD.TUMOR.FREE
table:
Table continues below X0.ALIVE.OR.DEAD.TUMOR.FREE X0.ALIVE.OR.DEAD.TUMOR.FREE 161 X1.DEAD.WITH.TUMOR 24 X1.DEAD.WITH.TUMOR X0.ALIVE.OR.DEAD.TUMOR.FREE 33 X1.DEAD.WITH.TUMOR 22 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.7625 0.2868 0.7035 0.8149 0.7708 AccuracyPValue McnemarPValue 0.6538 0.2893 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 0.8703 0.4 0.8299 0.4783 0.8299 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 0.8703 0.8496 0.7708 0.6708 0.8083 Balanced Accuracy 0.6351 mode: sens_spec
dots:
Matriz de Confusão:
- 161 previsões corretas para “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 22 previsões corretas para “X1.DEAD.WITH.TUMOR”.
- 33 erros ao classificar “X1.DEAD.WITH.TUMOR” como “XO.ALIVE.OR.DEAD.TUMOR.FREE”.
- 24 erros ao classificar “XO.ALIVE.OR.DEAD.TUMOR.FREE” como “X1.DEAD.WITH.TUMOR”
- Modelo tem boa capacidade de prever a classe “XO.ALIVE.OR.DEAD.TUMOR.FREE”, mas tem dificuldades em prever corretamente “X1.DEAD.WITH.TUMOR”.
- Embora treinado com Smote, o modelo pode estar enviesado para prever “XO.ALIVE.OR.DEAD.TUMOR.FREE” com mais frequência.
Métricas de Desempenho:
- Accuracy: 76.25% → Classifica corretamente na maioria dos casos.
- Kappa: 0.2886 → Indica um acordo moderado entre as previsões e os valores reais.
- Especificidade: 0.4 → Tem dificuldade em identificar corretamente os casos negativos.
- Balanced Accuracy: 0.6351 → O modelo não está completamente equilibrado, pois a especificidade é baixa.
# Previsão de probabilidades
pred_prob <- suppressWarnings(predict(fit_xgb1, testData, type = "prob"))
correct_col <- grep("ALIVE.OR.DEAD.TUMOR.FREE", colnames(pred_prob), value = TRUE)
f1_smote <- F1_Score(y_pred = pred_class, y_true = testData$SurvivalStatus, positive = "ALIVE.OR.DEAD.TUMOR.FREE")
# Precision, Recall
positive_class <- "X0.ALIVE.OR.DEAD.TUMOR.FREE"
precision <- posPredValue(pred_class, testData$SurvivalStatus, positive = positive_class)
recall <- sensitivity(pred_class, testData$SurvivalStatus, positive = positive_class)
f1 <- 2 * (precision * recall) / (precision + recall)
cat("Recall:", recall, "\n")Recall: 0.8702703
Precision: 0.8298969
F1 Score: 0.85
- Recall: 0.8703 → Identifica 87.03% dos casos positivos corretamente.
- Precisão: 0.8299 → Quando prevê “XO.ALIVE.OR.DEAD.TUMOR.FREE”, está correto 82.99% das vezes.
- F1 Score: 0.85 → Indica um bom equilíbrio entre precisão e recall.
# ROC
roc_curve <- roc(response = testData$SurvivalStatus, predictor = pred_prob[, correct_col])
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.7493857
Forma da Curva:
- A curva ROC está próxima do canto superior esquerdo, indicando alta sensibilidade e especificidade. Isso significa que o modelo classifica corretamente a maioria dos exemplos, minimizando falsos positivos e falsos negativos.
Comparação com um Classificador Aleatório
- A linha diagonal representa um classificador aleatório (AUC = 0.5).
- Como a curva ROC do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score
- AUC: 0.7494 → O modelo tem capacidade razoável de distinguir entre as classes, mas não é perfeito (abaixo de 0.8).
# Importância das features
importance <- varImp(fit_xgb1, scale = FALSE)
# Extrair variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
90187 |
90187 |
0.1019811 |
85236 |
85236 |
0.0785528 |
6427 |
6427 |
0.0456430 |
23590 |
23590 |
0.0373581 |
7060 |
7060 |
0.0370359 |
84236 |
84236 |
0.0332272 |
80311 |
80311 |
0.0293744 |
22944 |
22944 |
0.0273073 |
26959 |
26959 |
0.0257791 |
64072 |
64072 |
0.0245815 |
5116 |
5116 |
0.0226117 |
9770 |
9770 |
0.0222662 |
7764 |
7764 |
0.0203294 |
358 |
358 |
0.0199682 |
9862 |
9862 |
0.0191263 |
7289 |
7289 |
0.0187641 |
57653 |
57653 |
0.0183665 |
5021 |
5021 |
0.0177597 |
222161 |
222161 |
0.0140640 |
81575 |
81575 |
0.0122191 |
Importância dos Genes:
- Gene mais relevante: 90187, com um peso de 0.1019811. Este gene tem um impacto significativo na classificação do modelo.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.03), sugerindo que são os principais responsáveis pela decisão do modelo
# Boxplot das probabilidades
boxplot(pred_prob[, correct_col] ~ testData$SurvivalStatus,
col = palette_3,
main = "Class Probability Distribution", ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe “X0.ALIVE.OR.DEAD.TUMOR.FREE” (vermelho):
- A maioria das previsões tem probabilidades próximas de 1, indicando alta confiança do modelo nesta classe.
- Existem alguns outliers com probabilidades mais baixas (~0.2), sugerindo casos ambíguos.
Classe “X1.DEAD.WITH.TUMOR” (azul):
- A distribuição é mais espalhada, com valores variando de 0.2 a 0.8.
- A mediana está próxima de 0.5, indicando que o modelo não tem tanta certeza ao prever esta classe. Pode sugerir dificuldade na separação entre as classes, levando a previsões menos confiáveis.
## C. XGBoost com SMOTE - Sex
# Preparar dados
clean_data_for_2part <- clini_data_no_na
# Definir variável resposta
f_sex <- as.factor(clean_data_for_2part$sex)
names(f_sex) <- clean_data_for_2part$sample_ID
# Alinhar dados de expressão com amostras com dados clínicos
common_samples <- intersect(colnames(expr_all), names(f_sex))
logCPM_sub <- expr_all[, common_samples]
labels <- f_sex[common_samples]
# Construir dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub)) # genes como colunas
data_model$Sex <- factor(labels)
# Ajustar levels para uso com SMOTE e xgboost
levels(data_model$Sex) <- make.names(levels(data_model$Sex))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Sex, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# SMOTE
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
smote_result <- SMOTE(trainData_numeric, trainData$Sex, K = 3, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "Sex"
smote_data$Sex <- factor(smote_data$Sex, levels = levels(trainData$Sex))
# Treinar modelo XGBoost
set.seed(123)
capture.output({
fit_xgb2 <- train(Sex ~ ., data = smote_data, method = "xgbTree",
trControl = ctrl,
metric = "ROC", tuneLength=3)
}, file = "NUL")
# Avaliação no conjunto de teste
pred_class <- predict(fit_xgb2, testData)
conf_matrix <- confusionMatrix(pred_class, testData$Sex)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive: Female
table:
Female Male Female 106 3 Male 0 131 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.9875 0.9747 0.9639 0.9974 0.5583 AccuracyPValue McnemarPValue 2.056e-55 0.2482 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Neg Pred Value Precision 1 0.9776 0.9725 1 0.9725 Table continues below Recall F1 Prevalence Detection Rate Detection Prevalence 1 0.986 0.4417 0.4417 0.4542 Balanced Accuracy 0.9888 mode: sens_spec
dots:
Matriz de confusão:
- 106 previsões corretas para “Female” e 131 previsões corretas para “Male”.
- Apenas 3 erros ao classificar “Female” como “Male”.
- Nenhum erro ao classificar “Male”.
Modelo tem alta precisão e sensibilidade, especialmente para a classe “Female”.
- Treinado com SMOTE (Synthetic Minority Over-sampling Technique), que ajuda a lidar com desbalanceamento de classes.
Métricas de Desempenho:
- Accuracy: 98.75% → classifica corretamente na maioria dos casos.
- Kappa: 0.9747 → Excelente acordo entre as previsões e os valores reais.
- Especificidade: 0.9776 → Identifica quase todos os exemplos da classe “Male” bem.
# Previsão de probabilidades
pred_prob <- suppressWarnings(predict(fit_xgb2, testData, type = "prob"))
correct_col <- grep("Female", colnames(pred_prob), value = TRUE)
# Precision, Recall, F1
positive_class <- "Female"
precision <- posPredValue(pred_class, testData$Sex, positive = positive_class)
recall <- sensitivity(pred_class, testData$Sex, positive = positive_class)
f1 <- 2 * (precision * recall) / (precision + recall)
# Print precision, recall,F1
cat("Precision:", precision, "\n")Precision: 0.9724771
Recall: 1
F1 Score: 0.986
- Precisão: 0.9725 → Quando o modelo prevê “Female”, está correto 97.25% das vezes.
- Recall: 1 → Identifica todos os exemplos da classe “Female” corretamente.
- F1 Score: 0.986 → Indica um excelente equilíbrio entre precisão e recall.
# Curva ROC e AUC
roc_curve <- roc(response = testData$Sex, predictor = pred_prob[, correct_col])
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.9997888
Forma da Curva:
- A curva está muito próxima do canto superior esquerdo, indicando alta sensibilidade e especificidade, minimizando falsos positivos e falsos negativos.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score:
- AUC: 0.9998 → Quase perfeita separação entre as classes “Female” e “Male”. Distingue muito bem entre as duas classes.
# Importância das features
importance <- varImp(fit_xgb2, scale = FALSE)
# Extrair variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
22829 |
22829 |
0.7882625 |
8242 |
8242 |
0.0876079 |
7403 |
7403 |
0.0294961 |
6622 |
6622 |
0.0088017 |
79192 |
79192 |
0.0076825 |
1757 |
1757 |
0.0061126 |
51123 |
51123 |
0.0060821 |
1627 |
1627 |
0.0052512 |
388336 |
388336 |
0.0051415 |
285966 |
285966 |
0.0051263 |
1513 |
1513 |
0.0051139 |
94056 |
94056 |
0.0045483 |
146198 |
146198 |
0.0043606 |
11251 |
11251 |
0.0036224 |
25879 |
25879 |
0.0035475 |
692312 |
692312 |
0.0034098 |
399687 |
399687 |
0.0033947 |
57619 |
57619 |
0.0030025 |
8226 |
8226 |
0.0026597 |
23438 |
23438 |
0.0025830 |
Importância dos Genes:
- Gene mais relevante: 22829 com um peso de 0.788.
- Outros genes com menor impacto incluem 8242 (0.0876) e 7403 (0.0295).
- A maioria dos genes tem um impacto relativamente baixo, exceto o primeiro.
Isso sugere que alguns genes têm um papel muito forte na classificação, enquanto outros contribuem menos.
# Boxplot das probabilidades
boxplot(pred_prob[, correct_col] ~ testData$Sex,
col = palette_3,
main = "Distribuição das Probabilidades - Sexo", ylab = "Probabilidade Prevista")
Boxplot - Distribuição das Probabilidades:
Classe “Female” (vermelho):
- A maioria das previsões tem probabilidades próximas de 1, indicando alta confiança do modelo.
Classe “Male” (azul):
- A maioria das previsões tem probabilidades próximas de 0, também indicando alta confiança.
- Poucos outliers.
Isso confirma que o modelo classifica com alta confiança, mas pode haver casos ambíguos.
## C. XGBoost com SMOTE - Ancestralidade
# Preparar dados
clean_data_for_2part <- clini_data_no_na
# Preparar f_ancestralidade sem NA
f_ancestralidade <- as.factor(clean_data_for_2part$ancestry)
names(f_ancestralidade) <- clean_data_for_2part$sample_ID
f_ancestralidade <- f_ancestralidade[!is.na(f_ancestralidade)]
# Alinhar com dados de expressão
common_samples <- intersect(colnames(expr_all), names(f_ancestralidade))
logCPM_sub <- expr_all[, common_samples]
labels <- f_ancestralidade[common_samples]
# Criar dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub)) # genes como colunas
data_model$Ancestralidade <- factor(labels)
data_model <- na.omit(data_model) # Garantir que não há NAs
levels(data_model$Ancestralidade) <- make.names(levels(data_model$Ancestralidade))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$Ancestralidade, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# Remover classes com poucas amostras no treino
min_class_size <- 5
keep_classes <- names(which(table(trainData$Ancestralidade) >= min_class_size))
trainData <- trainData[trainData$Ancestralidade %in% keep_classes, ]
trainData$Ancestralidade <- droplevels(trainData$Ancestralidade)
# Aplicar SMOTE
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
smote_result <- SMOTE(trainData_numeric, trainData$Ancestralidade, K = 3, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "Ancestralidade"
smote_data$Ancestralidade <- factor(smote_data$Ancestralidade, levels = levels(trainData$Ancestralidade))
# Treinar modelo XGBoost
set.seed(123)
capture.output({
fit_xgb3 <- train(Ancestralidade ~ ., data = smote_data, method = "xgbTree",
trControl = ctrl1,
metric = "ROC",
tuneLength=3)
}, file = "NUL")
# Avaliação no conjunto de teste
pred_class <- predict(fit_xgb3, testData)
conf_matrix <- confusionMatrix(pred_class, testData$Ancestralidade)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive:
table:
AFR AFR_ADMIX AMR EAS EUR EUR_ADMIX SAS SAS_ADMIX AFR 0 0 0 0 0 0 0 0 AFR_ADMIX 0 0 0 0 0 0 0 0 AMR 0 0 0 0 0 0 0 0 EAS 0 0 0 0 0 0 0 0 EUR 7 4 2 4 218 3 0 1 EUR_ADMIX 0 0 0 0 0 0 0 0 SAS 0 0 0 0 0 0 0 0 SAS_ADMIX 0 0 0 0 0 0 0 0 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.9121 0 0.8688 0.9448 0.9121 AccuracyPValue McnemarPValue 0.5577 NA byClass:
Table continues below Sensitivity Specificity Pos Pred Value Class: AFR 0 1 NA Class: AFR_ADMIX 0 1 NA Class: AMR 0 1 NA Class: EAS 0 1 NA Class: EUR 1 0 0.9121 Class: EUR_ADMIX 0 1 NA Class: SAS NA 1 NA Class: SAS_ADMIX 0 1 NA Table continues below Neg Pred Value Precision Recall F1 Class: AFR 0.9707 NA 0 NA Class: AFR_ADMIX 0.9833 NA 0 NA Class: AMR 0.9916 NA 0 NA Class: EAS 0.9833 NA 0 NA Class: EUR NA 0.9121 1 0.954 Class: EUR_ADMIX 0.9874 NA 0 NA Class: SAS NA NA NA NA Class: SAS_ADMIX 0.9958 NA 0 NA Table continues below Prevalence Detection Rate Detection Prevalence Class: AFR 0.02929 0 0 Class: AFR_ADMIX 0.01674 0 0 Class: AMR 0.008368 0 0 Class: EAS 0.01674 0 0 Class: EUR 0.9121 0.9121 1 Class: EUR_ADMIX 0.01255 0 0 Class: SAS 0 0 0 Class: SAS_ADMIX 0.004184 0 0 Balanced Accuracy Class: AFR 0.5 Class: AFR_ADMIX 0.5 Class: AMR 0.5 Class: EAS 0.5 Class: EUR 0.5 Class: EUR_ADMIX 0.5 Class: SAS NA Class: SAS_ADMIX 0.5 mode: sens_spec
dots:
Matriz de Confusão:
- 218 previsões corretas para “EUR”, mas teve dificuldades nas outras classes.
- Confirma que o modelo está enviesado para prever EUR corretamente, mas falha nas outras classes
Métricas de Desempenho:
- Accuracy: 91.21% → Classifica corretamente na maioria dos casos.
- AccuracyLower: 0.8688 e AccuracyUpper: 0.9448 → accuracy real do modelo pode variar nesse intervalo que representa uma boa estabilidade no desempenho
- AccuracyNull: 0.9121 → o modelo pode estar apenas a prever a classe dominante.
- Kappa: 0 → indica que o modelo não está a prever melhor do que um classificador aleatório.
- Especificidade: 0.4 → Tem dificuldade em identificar corretamente os casos negativos.
- Balanced Accuracy: 0.6351 → O modelo não está completamente equilibrado, pois a especificidade é baixa.
- p-value da accuracy (0.5577) indica que não há diferença estatística significativa entre o modelo e um classificador aleatório.
- McNemarPValue está como NA, o que pode indicar que o teste não foi realizado ou não é aplicável neste caso.
# Previsão de probabilidades (para AUC e boxplot)
pred_prob <- suppressWarnings(predict(fit_xgb3, testData, type = "prob"))
# Filtrar apenas as classes "EUR" e "AFR"
testData <- testData[testData$Ancestralidade %in% c("EUR", "AFR"), ]
pred_prob <- predict(fit_xgb3, testData, type = "prob")
# Escolher uma classe como "positiva" arbitrária (por exemplo, a primeira)
positive_class <- colnames(pred_prob)[1]
# Calcular Precision e Recall
precision_vec <- diag(conf_matrix$table) / colSums(conf_matrix$table)
recall_vec <- diag(conf_matrix$table) / rowSums(conf_matrix$table)
# Calcular Macro F1-score manualmente (sem yardstick)
f1_per_class <- function(conf_mat) {
by_class <- conf_mat$byClass
if (is.matrix(by_class)) {
f1 <- 2 * by_class[, "Sensitivity"] * by_class[, "Precision"] /
(by_class[, "Sensitivity"] + by_class[, "Precision"])
f1[is.na(f1)] <- 0
return(mean(f1))
} else {
# Only two classes
f1 <- 2 * by_class["Sensitivity"] * by_class["Precision"] /
(by_class["Sensitivity"] + by_class["Precision"])
return(ifelse(is.na(f1), 0, f1))
}
}
f1_macro <- f1_per_class(conf_matrix)
cat("F1 Score:", round(f1_macro, 3), "\n")F1 Score: 0.119
Recall: NaN NaN NaN NaN 0.912 NaN NaN NaN
Precision: 0 0 0 0 1 0 NaN 0
- Precisão: 1 → Prevê apenas a classe “EUR” com alta confiança.
- Recall: 0.912 → Identifica bem apenas a classe “EUR”.
- F1 Score: 0.119 → Indica um score muito baixo, o que indica um grande desequilíbrio entre precisão e recall.
roc_curve <- roc(response = testData$Ancestralidade,
predictor = pred_prob[, positive_class],
levels = c("EUR", "AFR"))
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.827654
# cálculo de AUC para multiclasse
multiclass_auc <- multiclass.roc(response = testData$Ancestralidade,
predictor = as.matrix(pred_prob))
auc_score <- auc(multiclass_auc)
cat("Multiclass AUC Score:", round(auc_score, 3), "\n")Multiclass AUC Score: 0.828
Forma da Curva:
- A curva segue um padrão escalonado, indicando que o modelo está a tomar decisões em diferentes limiares de classificação, mas está acima da linha diagonal, o que significa que o modelo tem um desempenho melhor do que um classificador aleatório.
Comparação com um Classificador Aleatório:
- A linha diagonal representa um classificador aleatório (AUC = 0.5). Como a curva do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score:
- AUC: 0.827654 → Distingue bem algumas das classes.
# Importância das features
importance <- varImp(fit_xgb3, scale = FALSE)
# Extração variavel importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
317749 |
317749 |
0.1255606 |
283345 |
283345 |
0.1217968 |
10901 |
10901 |
0.0838071 |
93100 |
93100 |
0.0751637 |
8418 |
8418 |
0.0720582 |
10352 |
10352 |
0.0527334 |
324 |
324 |
0.0512933 |
57536 |
57536 |
0.0507113 |
2530 |
2530 |
0.0487659 |
9732 |
9732 |
0.0443780 |
221442 |
221442 |
0.0402516 |
283726 |
283726 |
0.0282638 |
9652 |
9652 |
0.0234975 |
165679 |
165679 |
0.0217786 |
80224 |
80224 |
0.0207446 |
6231 |
6231 |
0.0190774 |
23313 |
23313 |
0.0171899 |
285513 |
285513 |
0.0107812 |
57211 |
57211 |
0.0105151 |
83872 |
83872 |
0.0099357 |
Importância dos Genes:
- Genes mais relevantes: 317749 com um peso de 0.125560 e Gene 283345 com peso 0.1217968.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.07), sugerindo que são os principais responsáveis pela decisão do modelo.
- O modelo pode estar dependente de poucos genes-chave, enquanto outros têm impacto reduzido. Isso pode indicar que algumas características podem ser removidas sem afetar o desempenho.
# Boxplot das probabilidades para cada classe
positive_class <- colnames(pred_prob)[1]
boxplot(pred_prob[, positive_class] ~ testData$Ancestralidade,
col = palette_3,
main = paste("Probabilidade predita - Classe:", positive_class),
ylab = "Predicted Probability")
Boxplot - Distribuição das Probabilidades:
Classe AFR:
- A mediana da probabilidade predita está próxima de 0.05.
- O intervalo interquartil (IQR) varia entre aproximadamente 0.02 e 0.1, indicando baixa dispersão.
- Não há outliers, sugerindo que o modelo prevê esta classe de forma consistente.
Outras classes
- EUR (Europeia) apresenta vários outliers, com probabilidades chegando a 0.35.
- As classes AFR_ADMIX, AMR, EAS, EUR_ADMIX, SAS e SAS_ADMIX têm probabilidades próximas de zero, indicando que o modelo raramente prevê essas classes como AFR.
## C. XGBoost com SMOTE - Cancer
# Preparar dados
clean_data_for_2part <- clini_data_no_na
# Definir variável resposta
f_cancer <- as.factor(clean_data_for_2part$cancer_type)
names(f_cancer) <- clean_data_for_2part$sample_ID
# Alinhar dados de expressão com amostras com dados clínicos
common_samples <- intersect(colnames(expr_all), names(f_cancer))
logCPM_sub <- expr_all[, common_samples]
labels <- f_cancer[common_samples]
# Construir dataframe para modelagem
data_model <- as.data.frame(t(logCPM_sub)) # genes como colunas
data_model$TipoTumor <- factor(labels)
# Ajustar levels para uso com SMOTE e xgboost
levels(data_model$TipoTumor) <- make.names(levels(data_model$TipoTumor))
# Dividir treino e teste
set.seed(123)
trainIndex <- createDataPartition(data_model$TipoTumor, p = 0.5, list = FALSE)
trainData <- data_model[trainIndex, ]
testData <- data_model[-trainIndex, ]
# SMOTE
trainData_numeric <- trainData[, sapply(trainData, is.numeric)]
smote_result <- SMOTE(trainData_numeric, trainData$TipoTumor, K = 3, dup_size = 2)
smote_data <- smote_result$data
colnames(smote_data)[ncol(smote_data)] <- "TipoTumor"
smote_data$TipoTumor <- factor(smote_data$TipoTumor, levels = levels(trainData$TipoTumor))
# Treinar modelo XGBoost
set.seed(123)
capture.output({
fit_xgb4 <- train(TipoTumor ~ ., data = smote_data, method = "xgbTree",
trControl = ctrl1,
metric = "ROC", tuneLength=3)
}, file = "NUL")
# Avaliação no conjunto de teste
pred_class <- predict(fit_xgb4, testData)
conf_matrix <- confusionMatrix(pred_class, testData$TipoTumor)
pander::pander(conf_matrix, caption = "Confusion Matrix Statistics (per class)")positive:
table:
Astrocytoma Oligoastrocytoma Oligodendroglioma Astrocytoma 58 20 14 Oligoastrocytoma 29 24 18 Oligodendroglioma 4 18 54 overall:
Table continues below Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull 0.569 0.3489 0.5036 0.6327 0.3808 AccuracyPValue McnemarPValue 2.91e-09 0.06554 byClass:
Table continues below Sensitivity Specificity Pos Pred Value Class: Astrocytoma 0.6374 0.7703 0.6304 Class: Oligoastrocytoma 0.3871 0.7345 0.338 Class: Oligodendroglioma 0.6279 0.8562 0.7105 Table continues below Neg Pred Value Precision Recall F1 Class: Astrocytoma 0.7755 0.6304 0.6374 0.6339 Class: Oligoastrocytoma 0.7738 0.338 0.3871 0.3609 Class: Oligodendroglioma 0.8037 0.7105 0.6279 0.6667 Table continues below Prevalence Detection Rate Class: Astrocytoma 0.3808 0.2427 Class: Oligoastrocytoma 0.2594 0.1004 Class: Oligodendroglioma 0.3598 0.2259 Detection Prevalence Balanced Accuracy Class: Astrocytoma 0.3849 0.7038 Class: Oligoastrocytoma 0.2971 0.5608 Class: Oligodendroglioma 0.318 0.7421 mode: sens_spec
dots:
Matriz de Confusão:
- Astrocytoma: 58 corretos, 20 classificados como Oligoastrocytoma, 14 como Oligodendroglioma.
- Oligoastrocytoma: 24 corretos, 29 classificados como Astrocytoma, 18 como Oligodendroglioma.
- Oligodendroglioma: 54 corretos, 18 classificados como Oligoastrocytoma, 4 como Astrocytoma.Astrocytoma e Astrocytoma e Oligodendroglioma têm melhor desempenho, com sensibilidade acima de 62%. Oligoastrocytoma tem o pior desempenho, com sensibilidade de apenas 38.71%, indicando que o modelo tem dificuldades em identificar corretamente essa classe.
Métricas de Desempenho:
- Accuracy: 56.9% → O modelo classifica corretamente 56.9% dos exemplos, o que pode ser melhorado.
- Kappa: 0.3489 → Indica concordância moderada entre previsões e valores reais.
- Mcnemar’s Test P-Value: 0.06554 → P-valor relativamente alto, indicando que o modelo pode estar enviesado para certas classes, sem grande distinção entre elas.
Especificidade: O modelo distingue bem as classes, mas ainda há erros na previsão de Oligoastrocytoma.
- Astrocytoma: 77.03%
- Oligoastrocytoma: 73.45%
- Oligodendroglioma: 85.62%
Balanced Accuracy: Oligoastrocytoma tem a pior acurácia balanceada, indicando que o modelo tem dificuldades em prever essa classe corretamente.
- Astrocytoma: 70.38%
- Oligoastrocytoma: 56.08%
- Oligodendroglioma: 74.21%
# Previsão de probabilidades
pred_prob <- suppressWarnings(predict(fit_xgb4, testData, type = "prob"))
# Escolher uma classe como "positiva" arbitrária (por exemplo, a primeira)
positive_class <- levels(testData$TipoTumor)[1]
roc_curve <-roc(response = testData$TipoTumor, predictor = pred_prob[, positive_class])
auc_score <- auc(roc_curve)
cat("AUC Score:", auc_score, "\n")AUC Score: 0.718894
# Curva ROC
plot(roc_curve, col = palette_3, main = "ROC Curve for XGBoost (com SMOTE) - TipoTumor")
Forma da Curva:
- A curva segue um padrão escalonado, indicando que o modelo está a tomar decisões em diferentes limiares de classificação.
Comparação com um Classificador Aleatório:
A linha diagonal representa um classificador aleatório (AUC = 0.5).
Como a curva do modelo está bem acima dessa linha, isso confirma que o modelo tem um desempenho significativamente melhor do que um classificador aleatório.
AUC Score: 0.7648 → Tem uma boa capacidade de distinguir entre as classes. Acima de 0.75 indica um modelo útil.
# Calcular F1-score, Recall e Precision (para a mesma classe positiva)
f1_smote <- F1_Score(y_pred = pred_class, y_true = testData$TipoTumor, positive = positive_class)
precision_vec <- diag(conf_matrix$table) / colSums(conf_matrix$table)
recall_vec <- diag(conf_matrix$table) / rowSums(conf_matrix$table)
cat("F1 Score:", f1_smote, "\n")F1 Score: 0.6338798
Recall: 0.63 0.338 0.711
Precision: 0.637 0.387 0.628
Recall: O modelo identifica bem Astrocytoma e Oligodendroglioma, mas tem dificuldades com Oligoastrocytoma.
- Astrocytoma: 63.74%
- Oligoastrocytoma: 38.71%
- Oligodendroglioma: 62.79%
Precision:
Astrocytoma: 77.55%
Oligoastrocytoma: 77.38%
Oligodendroglioma: 80.37%
F1 Score: 0.6857 → O F1 Score mede o equilíbrio entre precisão e recall. Um valor de 0.6857 indica um desempenho razoável.
# Importância das features
importance <- varImp(fit_xgb4, scale = FALSE)
# Extrair variável importance como data frame
importance_df <- as.data.frame(importance$importance)
importance_df$Gene <- rownames(importance_df)
importance_df <- importance_df[order(-importance_df$Overall), ] # Sort by importance
# Seleção top 20
top_20 <- head(importance_df, 20)
# Tabela
knitr::kable(top_20[, c("Gene", "Overall")], caption = "Top 20 Most Important Genes")| Gene | Overall | |
|---|---|---|
55616 |
55616 |
0.0452976 |
1436 |
1436 |
0.0253459 |
9776 |
9776 |
0.0227275 |
57558 |
57558 |
0.0214301 |
8915 |
8915 |
0.0179006 |
1404 |
1404 |
0.0169052 |
285237 |
285237 |
0.0160138 |
1288 |
1288 |
0.0154437 |
23585 |
23585 |
0.0146489 |
58529 |
58529 |
0.0129332 |
3587 |
3587 |
0.0117539 |
90141 |
90141 |
0.0117313 |
2203 |
2203 |
0.0108208 |
10123 |
10123 |
0.0103094 |
5806 |
5806 |
0.0102897 |
219699 |
219699 |
0.0101287 |
56160 |
56160 |
0.0095977 |
5432 |
5432 |
0.0095320 |
219348 |
219348 |
0.0090669 |
494514 |
494514 |
0.0090427 |
Importância dos Genes:
- Genes mais relevantes: 55616 com um peso de 0.0452976 e Gene 1436 com peso 0.0253459.
- Os primeiros 5 genes têm valores relativamente altos (acima de 0.02), sugerindo que são os principais responsáveis pela decisão do modelo.
- O modelo pode estar dependente de poucos genes-chave, enquanto outros têm impacto reduzido. Isso pode indicar que algumas características podem ser removidas sem afetar o desempenho.
# Boxplot das probabilidades
boxplot(pred_prob[, positive_class] ~ testData$TipoTumor,
col = palette_3,
main = "Distribuição das Probabilidades - TipoTumor", ylab = "Probabilidade Prevista")
Boxplot - Distribuição das Probabilidades:
Astrocytoma:
- Mediana alta, indicando que o modelo prevê esta classe com confiança.
- Intervalo interquartil (IQR) relativamente estreito, sugerindo baixa variabilidade nas previsões.
- Poucos outliers, indicando que o modelo é consistente ao prever Astrocytoma.
Oligoastrocytoma:
- Mediana mais baixa, indicando que o modelo tem dificuldades em prever esta classe corretamente.
- Maior dispersão, sugerindo incerteza nas previsões.
- Outliers presentes, indicando que algumas previsões são altamente erráticas.
Oligodendroglioma:
- Mediana alta, semelhante à de Astrocytoma.
- Menos variabilidade do que Oligoastrocytoma, indicando que o modelo prevê esta classe com mais confiança.
- Poucos outliers, sugerindo boa estabilidade nas previsões.
Interpretação:
- O modelo tem alta confiança ao prever Astrocytoma e Oligodendroglioma, mas incerteza ao prever Oligoastrocytoma.
- A maior variabilidade nas previsões de Oligoastrocytoma pode estar impactando negativamente a acurácia geral do modelo.
12. Comparação dos modelos por variável
# Variável Survival
min_resamples <- min(sapply(list(fit_rf1, fit_svm1, fit_xgb1), function(x) nrow(x$resample)))
fit_rf1$resample <- fit_rf1$resample[order(fit_rf1$resample$Resample), ]
fit_svm1$resample <- fit_svm1$resample[order(fit_svm1$resample$Resample), ]
fit_xgb1$resample <- fit_xgb1$resample[order(fit_xgb1$resample$Resample), ]
# Comparar resultados
results_survival <- resamples(list(RF = fit_rf1, SVM = fit_svm1, XGB = fit_xgb1))
summary_df <- as.data.frame(summary(results_survival)$statistics)
summary_df <- rownames_to_column(summary_df, "Metric_Model")
datatable(summary_df,
options = list(pageLength = 10, autoWidth = TRUE),
caption = "Summary of Resampling Metrics")
Análise da comparação:
Accuracy:
SVM tem a maior mediana de acurácia (0.9429), indicando que, em média, este é o modelo mais preciso.
XGB e RF têm medianas semelhantes (0.9143), mas RF tem um mínimo mais baixo (0.8571), sugerindo maior variabilidade.
SVM e XGB atingem a mesma accuracy máxima (0.9857), indicando que ambos podem ter um desempenho bom em alguns casos.
Kappa Score:
SVM tem a maior mediana de Kappa (0.8861), indicando que ele tem a melhor concordância entre previsões e valores reais.
XGB tem um Kappa semelhante ao RF na mediana (0.8291 vs. 0.8280), mas atinge um máximo igual ao SVM (0.8863).
RF tem um Kappa mínimo mais baixo (0.7143), sugerindo que pode ter mais variação na concordância.
# Variável Sex
min_resamples <- min(sapply(list(fit_rf2, fit_svm2, fit_xgb2), function(x) nrow(x$resample)))
fit_rf2$resample <- fit_rf2$resample[1:min_resamples, ]
fit_svm2$resample <- fit_svm2$resample[1:min_resamples, ]
fit_xgb2$resample <- fit_xgb2$resample[1:min_resamples, ]
# Comparar resultados
results_sex <- resamples(list(RF = fit_rf2, SVM = fit_svm2, XGB = fit_xgb2))
summary_df_sex <- as.data.frame(summary(results_sex)$statistics)
summary_df_sex <- rownames_to_column(summary_df_sex, "Metric_Model")
datatable(summary_df_sex,
options = list(pageLength = 10, autoWidth = TRUE),
caption = "Summary of Resampling Metrics (Sex)")
Análise de comparação:
Accuracy:
XGB: Maior acurácia entre os três modelos, com valores próximos de 1. Baixa variabilidade, indicando alta estabilidade nas previsões. O ponto médio (verde) está próximo do limite superior, sugerindo alta consistência.
SVM: Acurácia moderada, em torno de 0.8. Maior variabilidade do que XGB, indicando flutuações no desempenho. O ponto médio (azul) está abaixo do limite superior, sugerindo alguma inconsistência.
RF: Menor acurácia entre os três modelos, em torno de 0.7. Maior variabilidade, indicando instabilidade nas previsões. O ponto médio (vermelho) está no centro do boxplot, sugerindo desempenho inconsistente.
Kappa:
XGB tem o maior Kappa, com valores próximos de 1, indicando alta concordância entre previsões e valores reais.
SVM tem um Kappa moderado, com valores em torno de 0.5, sugerindo uma concordância razoável, mas inferior ao XGB.
RF tem o menor Kappa, com valores próximos de 0.3, indicando baixa concordância entre previsões e valores reais.
# Variável Ancestralidade
min_resamples <- min(sapply(list(fit_rf3, fit_svm3, fit_xgb3), function(x) nrow(x$resample)))
fit_rf3$resample <- fit_rf3$resample[1:min_resamples, ]
fit_svm3$resample <- fit_svm3$resample[1:min_resamples, ]
fit_xgb3$resample <- fit_xgb3$resample[1:min_resamples, ]
# Comparar resultados
results_ancestralidade <- resamples(list(RF = fit_rf3, SVM = fit_svm3, XGB = fit_xgb3))
summary_df_ancestralidade <- as.data.frame(summary(results_ancestralidade)$statistics)
summary_df_ancestralidade <- rownames_to_column(summary_df_ancestralidade, "Metric_Model")
datatable(summary_df_ancestralidade,
options = list(pageLength = 10, autoWidth = TRUE),
caption = "Summary of Resampling Metrics (Ancestralidade)")# Comparação Visual
bwplot(results_ancestralidade,
col = palette_3,
main = "Medidas ", ylab = "Model")
Comparação Geral dos Modelos:
XGB (XGBoost): Melhor desempenho na maioria das métricas, especialmente AUC, Balanced Accuracy e F1 Score. Menor variabilidade, indicando maior estabilidade nas previsões. Alta precisão e recall, sugerindo boa separação entre classes.
SVM (Support Vector Machine): Desempenho intermediário, com boa precisão e recall, mas maior variabilidade. Pode ser menos estável do que XGB, mas ainda apresenta bons resultados. Kappa e Balanced Accuracy são razoáveis, indicando concordância moderada.
RF (Random Forest): Menor desempenho em várias métricas, especialmente Kappa e logLoss. Maior variabilidade, sugerindo instabilidade nas previsões. Pode ser útil para interpretação, mas não é o melhor modelo em termos de precisão.
Análise das Métricas:
AUC e Balanced Accuracy: XGB tem os melhores valores, indicando excelente separação entre classes. SVM tem desempenho razoável, mas com maior variabilidade. RF tem os piores valores, sugerindo dificuldade na separação das classes.
F1 Score e Recall: XGB tem o melhor F1 Score, indicando bom equilíbrio entre precisão e recall. SVM tem recall alto, mas pode ter menor precisão. RF tem menor F1 Score, sugerindo baixa concordância entre previsões e valores reais.
Kappa e logLoss: XGB tem o melhor Kappa, indicando alta concordância entre previsões e valores reais. SVM tem um Kappa razoável, mas com maior variabilidade. RF tem o menor Kappa e maior logLoss, sugerindo baixa confiabilidade.
# Variável Cancer
min_resamples <- min(sapply(list(fit_rf4, fit_svm4, fit_xgb4), function(x) nrow(x$resample)))
fit_rf4$resample <- fit_rf4$resample[1:min_resamples, ]
fit_svm4$resample <- fit_svm4$resample[1:min_resamples, ]
fit_xgb4$resample <- fit_xgb4$resample[1:min_resamples, ]
# Comparar resultados
results_cancer <- resamples(list(RF = fit_rf4, SVM = fit_svm4, XGB = fit_xgb4))
summary_df_cancer <- as.data.frame(summary(results_cancer)$statistics)
summary_df_cancer <- rownames_to_column(summary_df_cancer, "Metric_Model")
datatable(summary_df_cancer,
options = list(pageLength = 10, autoWidth = TRUE),
caption = "Summary of Resampling Metrics (Cancer)")
Comparação Geral dos Modelos:
XGB (XGBoost): Melhor desempenho na maioria das métricas, especialmente AUC, Balanced Accuracy e F1 Score. Menor variabilidade, indicando maior estabilidade nas previsões. Alta precisão e recall, sugerindo boa separação entre classes.
SVM (Support Vector Machine): Desempenho intermediário, com boa precisão e recall, mas maior variabilidade. Pode ser menos estável do que XGB, mas ainda apresenta bons resultados Kappa e Balanced Accuracy são razoáveis, indicando concordância moderada.
RF (Random Forest): Menor desempenho em várias métricas, especialmente Kappa e logLoss. Maior variabilidade, sugerindo instabilidade nas previsões. Pode ser útil para interpretação, mas não é o melhor modelo em termos de precisão.
Análise das Métricas:
AUC e Balanced Accuracy: XGB tem os melhores valores, indicando excelente separação entre classes. SVM tem desempenho razoável, mas com maior variabilidade. RF tem os piores valores, sugerindo dificuldade na separação das classes.
F1 Score e Recall: XGB tem o melhor F1 Score, indicando bom equilíbrio entre precisão e recall. SVM tem recall alto, mas pode ter menor precisão. RF tem menor F1 Score, sugerindo baixa concordância entre previsões e valores reais.
Kappa e logLoss: XGB tem o melhor Kappa, indicando alta concordância entre previsões e valores reais. SVM tem um Kappa razoável, mas com maior variabilidade. RF tem o menor Kappa e maior logLoss, sugerindo baixa confiabilidade.